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ABSTRACT 


This  research  seeks  to  develop  non-invasive  burn  depth  evaluation  methods  from 
non-contacting  visible  and  near-infrared  spectroscopic  measurements.  In  previous  years, 
we  demonstrated  that  features  of  the  optical  reflection  spectra  of  burn  wounds  can  be 
correlated  with  the  depth  of  burn.  An  imaging  system  was  built  which  determined,  with 
accuracy  equal  to  or  better  than  that  of  a  skilled  burn  surgeon,  the  probability  that  burn 
sites  would  heal  within  three  weeks  from  date  of  Injury.  Our  goal  for  the  current  project 
is  to  investigate  the  optical  reflectance  properties  of  burns,  utilizing  the  techniques  of 
multivariate  analysis,  in  order  to  improve  the  reliability  of  this  instrument. 

Excellent  progress  has  been  made  toward  achievement  of  the  goals  of  this  project. 
In  the  first '  ar  a  commercial  spectrophotometer  was  purchased  and  modified  to  make 
optical  reflection  measurements  in  both  the  visible  and  near-infrared  regions  (450-1800 
nm).  A  library  of  reference  spectra  was  acquired  with  this  instrument. 

During  the  second  year,  experiments  were  conducted  to  understand  the  major 
components  of  the  reflectance  spectra  of  human  skin  /n  vivo.  Two  dynamic  processes 
which  influence  in  vivo  spectra  were  studied  to  improve  our  knowledge  of  burn 
physiology:  temperature  and  ischemia.  These  were  modeled  in  healthy  human  subjects. 
The  spectrum  of  pure  water,  which  is  highly  temperature  dependent  and  dominates  the 
spectra  of  biological  tissues,  was  mathematically  analyzed  for  the  purpose  of 
compensating  for  water  absorbances  in  tissues. 

Thirty  four  patients  at  Harborview  Burn  Center  were  studied,  and  the  reflectance 
spectra  of  their  burns  analyzed  with  multivariate  statistics.  An  unexpected  absorption 
band  correlating  to  burn  depth  was  observed  at  630  nm  and  was  identified  as  arising  from 
methemoglobin.  Based  on  the  foregoing  results,  a  more  reliable  algorithm  was  developed 
for  the  imaging  system.  The  new  algorithm  incorporates  wavelengths  which  measure  this 
methemoglobin  absorption.  It  has  been  tested  by  simulating  imaging  system  responses 
from  spectrophotometer  data.  The  new  model  predicts  the  healing  potential  of  110  burn 
sites  with  88%  accuracy,  significantly  better  than  simulated  responses  from  the  old 
imaging  system  (79%). 
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INTRODUCTION 

Burn  injuries  affect  more  than  two  million  American  victims  each  year.  Two 
hundred  thousand  of  these  will  require  hospitalization  and  10,000  will  die.’  Skin  grafting 
has  significantly  reduced  infection,  mortality,  and  disfigurement  in  patients  with  deep 
wounds.  Shallow  wounds,  however,  produce  smoother  skin  when  they  are  allowed  to 
heal  on  their  own.  As  a  general  rule,  only  burns  which  do  not  heal  in  21  days  benefit 
from  grafting^’^.  An  accurate  diagnosis  may  be  obtained  by  waiting  for  3  weeks  before 
deciding  to  operate,  but  it  has  been  shown  that  grafting  within  the  first  few  days  produces 
better  surgical  results  and  reduces  hospital  stay*^.  This  has  been  the  driving  force  behind 
the  development  of  methods  to  predict  potential  healing  times. 

For  a  wound  to  heal  well  (and  within  21  days),  it  must  re-epithelialize  from 
remnants  of  epithelium  in  the  lining  of  ducts  and  hair  follicles.  No  one  has  yet  proposed 
a  method  to  predict  healing  by  measuring  the  abundance  of  epithelium.  Instead,  most 
estimate  the  depth  of  irreparably  damaged  tissue  with  respect  to  the  total  thickness  of  the 
dermis. 

Many  hospitals  today  use  a  technique  called  clinical  assessment  which  combines 
visual  inspection  with  simple  tests,  such  as  the  wound’s  sensitivity  to  touch'’.  Data 
collection,  even  for  the  largest  and  most  variegated  wounds,  requires  less  than  15 
minutes  and  does  not  cause  excessive  pain,  neither  from  the  tests  themselves  or  from 
prolonged  exposure  of  the  wound.  However,  the  accuracy  of  this  technique  is  highly 
dependent  on  the  experience  of  the  assessor.  Other  techniques,  including 
histopathological  inspection,®®’^  topical  application  of  dyes,®  radio  isotope  studies,® 
ultrasound,’®  thermography,”’’^pulse  oximetry,’®  laser  doppler  flowmetry,’^’® 
and  multispectral  photographic  analyzers,’®  ’^  have  proven  to  be  slower,  more  painful, 
or  less  reliable  and  are  less  commonly  used. 
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The  Bum  Depth  Indicator 

In  1980  we  built  the  Burn  Depth  Indicator  (BDI),  a  hand-held  probe  which 
illuminated  a  3  cm^  area  of  wound  with  pulsed  light  from  three  LED  chips  and  measured 
the  reflected  light.  This  technique  proved  to  be  as  fast  as  clinical  assessment,  less 
painful,  and  significantly  more  reliable.’®'^®  We  developed  a  similar  instrument  in  1986, 
which  provided  spatially  resolved  burn  wound  diagnosis.  This  device,  the  Imaging  Burn 
Depth  Indicator  (IBDI),  produces  a  color-coded  Image  of  a  wound.^'^’’^ 

Selection  of  the  wavelengths  used  by  the  BDI  (green  (550  nm),  red  (640  nm),  and 
near  infrared  (880  nm))  was  based  on  the  visible  appearance  of  burn  wounds  and  the 
discovery  that  near  infrared  light  is  reflected  more  by  deep  burns  than  shallow.*’®  Van 
Uew  interpreted  the  BDI  measurements  in  terms  of  blood  volume,  blood  oxygenation,  and 
eschar  thickness.  The  results  of  his  model  agreed  well  with  both  the  measured  data  and 
the  expected  physiology.*’®’^ 

Spectroscopic  Oximetry  in  vivo 

Measurement  of  blood  oxygenation  by  spectral  analysis  of  reflected  or  transmitted 
light  is  a  well  established  method.  Visible  (400-700  nm)  and  short-wave  near  infrared 
(700-1100  nm)  light  have  been  employed.  Hemoglobins  have  tt  tt*,  d  -  d*  and  f  -♦  f* 
electronic  transitions  in  this  region.  Oxidation  and  ligand  binding  influence  the  spectrum 
substantially  (Figure  1).  A  non-invasive  oximeter,  the  Pulsed  Oximeter,^'’  now  common 
in  hospitals,  uses  transmission  of  660  and  920  nm  light  through  a  finger  tip.  Its  results 
correlate  well  with  arterial  blood  Oj  saturation  measured  by  conventional,  invasive 
methods.®®  The  success  of  this  method  has  stimulated  research  on  other  non-invasive, 
spectrophotometric  oximeters.  These  new  instruments  measure  either  transmitted  or 
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reflected  light  and  have  been  applied  to  the  head,^®’^’^  brain, heart, and 
skin.*^®-^  These  methods  are  rarely  calibrated  to  give  absolute  values,  due  to  the 
unavailability  of  reference  data,  but  they  are  quite  useful  for  clinical  trend  monitoring  and 
kinetic  studies,  where  relative  changes  are  important. 

Near  Infrared  Specfroscopy 

Near  infrared  light  (700-2500  nm)  is  absorbed  by  other  biological  components 
besides  hemoglobins.  Most  near  infrared  (NIR)  absorbances  are  due  to  overtones  and 
combination  bands  of  vibrational  transitions,  primarily  C-H,  N-H,  and  0-H.  This  region 
has  been  used  in  the  agricultural  and  food  processing  industries  to  determine,  for 
instance,  protein,  oil,  and  moisture  content.®®®’  Samples  are  frequently  analyzed  in 
solid  or  powdered  form  by  diffuse  reflectance.  The  resulting  absorption  bands,  despite 
the  high  degree  of  light  scatter  in  the  samples,  are  often  linear  with  concentration. 
However  the  spectral  bands  are  broad  and  overlapping,  making  analysis  difficult. 
Calibrations  are  therefore  based  on  empirically  observed  spectral  behavior  rather  than  on 
first  principle  considerations. 

MultivariatB  Statistics 

The  key  that  has  unlocked  the  potential  of  NIR  spectroscopy  is  multivariate 
statistics,  mathematical  methods  for  analyzing  data  sets  which  include  multiple 
measurements  for  each  sample.  For  spectroscopy,  each  filter  or  each  monochromator 
setting  is  used  to  obtain  different  information,  producing  a  response  vector  for  each 
sample.  Multivariate  statistics  assists  the  analyst  to  find  an  empirical  mathematical  model 
which  will  best  reproduce  the  values  of  one  or  more  constituents,  given  the  response 
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vectors  of  a  set  of  calibration  samples.  Once  the  model  is  obtained,  it  may  be  used  to 
estimate  values  for  prediction  samples,  those  whose  constituents  are  unknown. 

Constituents  are  usually  chemical  concentrations,  but  may  also  be  physical 
properties,  e.g.  crystallinity  of  nylon  yarns,  or  abstract  properties,  e.g.  the  loaf  score  of 
flour  (the  volume,  appearance,  texture,  and  spring  of  bread  made  with  that  flour)  (31). 
The  property  needs  only  to  have  some  systematic  relation  to  the  spectra.  Multivariate 
statistics  can  be  used  to  analyze  the  spectra  of  burn  injuries  to  estimate  healing  potential, 
blood  oxygenation,  and  other  properties.  Further,  the  results  of  these  analyses  and  a 
more  precise  understanding  of  the  spectra  of  burned  skin  can  be  used  to  optimize  the 
algorithm  currently  used  by  the  IBDI. 
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PREUMINARY  RESEARCH 
THE  SPECTROMETER 

The  instrument  used  is  an  LT  Quantum  1200  spectrometer,  which  was  specifically 
modified  for  this  project  by  the  manufacturer.  It  employs  a  concave  holographic  grating 
and  a  manually  switched  pair  of  filters  to  scan  900-1800  nm  (first  order)  and  450-900  nm 
(second  order).  Light  is  conducted  to  the  sample  through  a  2  m  fiber  optic  cable  (0.5  cm 
bundle).  The  reflected  light  is  measured  by  two  45“  PbS  detectors  in  a  module  at  the 
end  of  the  cable  (Rgure  2).  The  module  is  suspended  1  to  2  cm  over  the  sample  and 
supported  by  a  tripod.  Data  is  collected  for  40  seconds  (100  scans)  and  analyzed  with 
software  written  in  the  MATLAB^  program  environment.  This  instrument  can  measure 
the  reflectance  of  1  cm^  of  skin  in  less  than  3  minutes,  including  the  time  required  to 
position  the  patient  and  detector.  It  is  completely  painless,  non-contacting,  and  non¬ 
destructive. 

REFERENCE  SPECTRA 

Reference  spectra  of  distilled  water,  mineral  oil,  collagen,  and  whole  blood  (Figures 
3-6)  were  obtained  using  the  LT  spectrometer  in  trans-reflectance  mode.  Band 
assignments  were  made  from  Colthup  type  charts  and  other  data*^  and  are  listed  in 
Table  I.  Reflectance  spectra  were  also  acquired  from  pork  fat,  pork  muscle,  beef  tendon, 
and  human  skin  in  vivo  (Figures  7-10).  These  spectra  are  composed  of  combinations  of 
hydrocarbon,  protein,  and  water  absorbances.  The  primary  absorber  in  biological  tissues 
is  water  (80%  by  weight).  Since  the  spectrum  of  water  is  highly  temperature  dependent, 
water  spectra  were  studied  in  detail  (pages  35"  ). 

Although  muscle,  tendon,  and  skin  are  primarily  water,  their  spectra  are  remarkably 
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different.  The  Beer-Lambert  law  (Absorbances  of  components  add  linearly), 

Aj)  =  -logiQ(/|j  /  /(j  jj)  =  bj  ()(£aj|(C)|,)  (1) 

Ajj  =  /“’  sample  absorbance  at  the  wavelength 
/,  jj  =  intensity  of  transmitted  radiation 
/q  ij  =  intensity  of  incident  radiation 

aj,,  =  absorptivity  of  the  /f“’  chemical  component  at  the  /'**’  wavelength, 
bj  =  path  length  of  light  through  sample  =  thickness  of  sample 
Cji,  =  concentration  of  the  /r”*  analyte  in  the  sample 

and  the  diffuse  reflectance  approximation 

A,j«log(1//?ij)  (2) 

«i|  =  (/r^.d.ii)/(W  (3) 

are  somewhat  inappropriate  for  these  samples;  the  assumptions  on  which  these 

equations  are  based  are  not  valid.  Light  scatter,  stray  light,  specular  reflectance,  and 
distribution  error  cause  pronounced  non-linearities  in  the  spectra  of  these  tissues.  The 
multiple  layers  of  tissues  in  skin  produce  additional  distortions.  These  non-linearities 
affect  both  the  visible  and  NIR  spectra  and  make  theoretical  analysis  extremely  difficult. 
On  the  other  hand,  they  may  provide  clues  to  physical  parameters,  such  as  surface 
texture  and  hydration,  distribution  of  blood  ceils  in  vessels  and  tissue,  disruption  of 
connective  tissue  layers,  etc.  These  parameters  generally  cannot  be  measured  directly, 
but  their  effects  on  the  spectra  of  burn  wounds  may  correlate  empirically  with  healing 
time. 
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IN  VIVO  STUDIES 
Skin  Temperature 

Preliminary  in  vivo  studies  were  done  to  observe  spectral  changes  due  to  expected 
dynamic  processes.  The  first  study  investigated  the  effect  of  temperature.  When  skin  is 
warmed,  the  capillaries  dilate  in  a  manner  similar  to  the  hyperemic  response  in  burns. 
Reflectance  spectroscopy  should  detect  temperature  dependent  shifts  in  water  bands  and 
an  increase  in  hemoglobin  concentration. 

Spectra  were  taken  from  the  forearms  of  seven  people  at  temperatures  ranging 
from  9  to  35  “C.  The  skin  was  cooled  with  water  and  heated  with  a  lamp.  Temperature 
was  measured  with  a  small  thermistor  at  the  surface  of  the  skin. 

One  subject  was  studied  in  detail  from  9  to  36  “C.  Figure  11a  shows  second 
derivative  NIR  spectra.  The  largest  variations  appear  at  1150  and  1213  nm.  Difference 
spectra  were  computed  using  the  spectrum  taken  at  9  “C  as  a  reference  (Figure  11b). 
The  large  peaks  at  970, 1155,  and  1390  nm  in  the  difference  spectra  are  attributed  to  the 
shifting  of  the  water  bands  at  elevated  temperatures  (Figure  11c). 

Similar  results  were  observed  in  all  seven  subjects.  The  subjects  were  divided  into 
two  groups.  Step-wise  linear  regression  was  used  to  build  two  3  wavelength  models 
using  each  group  as  a  separate  calibration  set.  Both  models  used  similar  wavelengths 
(1016,1073,1156  and  1001,1076,1150  nm)  The  1150,1156  and  1001,1016  nm 
wavelengths  correspond  to  spectral  changes  of  pure  water.  Each  model  was  used  to 
predict  temperatures  of  the  subjects  in  the  other  group  (Figure  12).  The  standard  errors 
of  prediction  (SEP)  were  1.3®  and  1.4®,  approximately  equal  to  the  estimated  accuracy 
of  the  reference  method  (1.5®C).  The  correlation  coefficients  were  both  0.982.  These 
results  suggest  that  it  may  be  possible  to  measure  the  surface  temperature  of  a  burn 
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wound  and  to  monitor  the  inflammation  of  the  wound  as  it  responds  to  injury  or  infection. 

Ischemia 

The  second  study  investigated  the  effect  of  ischemia,  which  may  cause  destruction 
in  the  zone  of  stasis.  The  largest  changes  would  be  expected  in  the  visible  region  where 
the  double-banded  absorbance  of  oxyhemoglobin  (Hb02)  would  be  replaced  by  the 
single,  broad  band  of  deoxyhemoglObin  (Hb)  (Figure  la). 

Ischemia  was  induced  in  the  forearm  of  a  healthy  volunteer  by  applying  a  pressure 
cuff  (150  mm  Hg)  for  8  minutes.  Spectra  were  recorded  in  the  visible  region  (450-780 
nm)  every  20  seconds,  while  circulation  was  restricted  and  for  8  minutes  after  it  was 
released. 

Rgure  13  shows  spectral  (second  derivative)  shifts  that  occurred  during  ischemia. 
The  two  bands  of  oxyhemoglobin  (540,  575  nm)  are  visible,  as  are  the  555  and  760  nm 
bands  of  deoxyhemoglobin.  Principal  Component  Analysis  (PCA)  resolves  the  spectra 
into  two  components  (Figures  14a,  14c)  corresponding  to  the  average  spectrum  (primarily 
HbOg)  and  the  Hb-HbOj  difference.  The  behavior  of  these  components  agrees  with  the 
known  physiology  of  cuff  induced  ischemia:  when  the  circulation  is  restricted,  the  blood 
becomes  deoxygenated.  When  the  restriction  is  released,  the  initial  surge  of  new  blood 
causes  a  hyperemic  increase  in  blood  volume,  primarily  oxygenated  arterial  blood,  which 
gradually  returns  to  normal  equilibrium.^  Figures  14b  and  14d  show  the  scores 
calculated  by  PCA,  the  estimated  "concentrations"  of  HbOj  and  Hb  components  that  were 
derived  from  the  spectra.  The  deoxygenation  (0-8  min.),  overshoot  (8  min.),  and 
reestablishment  of  equilibrium  (8-16  min.)  are  clearly  visible.  This  study  demonstrates  that 
the  LT  spectrometer  is  capable  of  detecting  changes  in  local  blood  oxygenation  in  the 
skin,  without  touching  the  patient  or  requiring  pulsed  blood  flow. 
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HUMAN  BURN  PATIENTS  AT  THE  NORTHWEST  BURN  CENTER 

The  main  body  of  this  research  project  is  the  analysis  and  interpretation  of  the 
spectra  acquired  from  the  skin  of  seriously  burned  patients  at  the  Northwest  Burn  Center 
at  Harborview  Hospital. 

EXPERIMENTAL 
Human  Subjects 

Forty-three  patients  were  studied  at  the  Northwest  Burn  Center,  Harborview 
Medical  Center  in  Seattle,  Washington,  over  a  13  month  period.  One  fifth  of  the  injuries 
were  selected  by  the  research  nurse  to  give  a  representative  sample  of  typical  superficial, 
shallow  partial,  deep  partial,  and  full  thickness  burns.  The  remainder  were  chosen 
because  they  were  indeterminant  or  unusual.  The  medically  fragile,  including  elderly, 
hemophiliac,  and  mentally  ill  patients,  were  excluded.  Children  were  also  excluded. 

Patients  were  examined  on  the  third  day  post-burn  as  was  done  in  previous 
studies.’®  A  few  of  the  unusual  injuries  were  studied  on  the  fourth  day,  because  the 
patients  were  transferred  from  other  hospitals  and  were  not  available  on  the  third  day. 

Two  to  five  sites  were  chosen  by  the  research  nurse,  often  with  the 
recommendation  of  the  attending  physician.  After  routine  cleaning  and  debridement 
(removal  of  loose  dead  tissue),  these  sites  were  wrapped  in  cellophane  to  prevent  drying 
and  infection.  Large  wounds  were  rebandaged,  except  for  the  selected  sites.  The  patient 
was  then  returned  to  his  own  roorh  and  positioned  either  in  bed  or  in  a  chair.  The 
cellophane  was  removed  prior  to  spectral  data  acquisition. 

One  spectrum  was  acquired  from  each  site,  using  the  LT  Quantum  1200 
spectrophotometer.®^  Site  locations  were  recorded  on  polaroid  photographs.  This 
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procedure  generally  took  15  to  30  minutes,  including  the  time  required  to  position  the 
patient  and  to  remove  the  equipment  afterwards.  If  the  patient  was  willing,  IBDI  images 
were  also  acquired,  generally  requiring  another  10  to  20  minutes. 

The  injuries  were  followed  by  visual  inspection  for  one  month.  The  sites  from 
which  spectra  had  been  acquired  were  classified  according  to  outcome  into  three 
categories:  Shallow  (healed  in  less  than  21  days).  Deep  (healed  in  more  than  21  days), 
and  Unknown.  Most  sites  which  were  grafted  were  classified  as  deep.  Those  which  were 
clearly  shallow  but  were  excised  due  to  their  proximity  to  a  more  severe  injury  were 
classified  as  shallow.  If  an  injury  was  not  grafted,  it  was  inspected  by  the  research  nurse 
on  or  near  the  21®*  day  to  determine  if  it  had  healed.  When  that  was  not  possible,  the 
information  was  obtained  by  phoning  the  patient.  Sites  for  which  there  was  insufficient 
information  were  classified  as  unknown  and  were  excluded  from  analysis. 

Of  the  sites  whose  outcomes  were  known,  34  were  deep  and  76  were  shallow. 
The  most  common  etiology  was  flame  (52%),  followed  by  scald  (22%),  chemical  (8%), 
flash  (8%),  contact  (6%),  and  grease  (4%).  The  mean  age  of  the  patients  was  33  years 
and  the  mean  TBSA  was  6%.  There  were  3  hispanics,  2  blacks,  and  1  oriental.  Burn 
locations  included  face,  arms,  hands,  palms,  backs,  sides,  legs,  and  feet. 

Mathematical  Methods 

The  multivariate  mathematical  methods  needed  to  interpret  data  include 
preprocessing,  calibration,  and  classification  techniques.  The  goal  is  to  find  a  useful 
relationship  between  a  sample’s  spectrum  (r)  and  the  values  of  its  constituents  (c). 
Constituents  are  usually  chemical  concentrations,  but  may  also  be  physical  properties, 
abstract  properties,  or  classifications.  The  property  needs  only  to  have  some  systematic 
relation  to  the  spectra. 
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The  notation  used  will  be  as  follows:  Matrices  will  be  written  in  bold,  upper  case 
(A,  B,  C),  row  vectors  in  bold,  lower  case  (a,  b,  c),  and  scalars  in  normal  italic  (A,  B,  C, 
a,  b,  c).  Column  vectors  will  be  written  as  the  transpose  of  row  vectors  (a^,  b^,  c^).  Each 
spectrum  will  be  referred  to  as  a  response  row  vector  r,  with  length  J,  where  J  is  the 
number  of  discrete  wavelengths  contained  in  each  spectrum.  /  spectra  grouped  into  a 
matrix  will  be  R,  with  dimensions  /  by  J.  The  elements  of  R  will  be  written  as  ry.  The 
matrix  C,  will  be  /  by  K,  where  K  is  the  number  of  spectrally  distinguishable  constituents. 
The  elements  of  C  are  Cj,^. 

Preprocessing 

Preprocessing  transforms  the  spectra  such  that  the  spectral  features  which 
correlate  to  the  constituents  are  maximized  and  the  remainder  are  minimized.  The  first 
step  is  usually  linearization  with  respect  to  constituent  concentration.  Normally  the 
logarithm  of  instrument  response  (R)  is  used  for  linearization: 

f?  =  /  /  /o  (3) 

where  Iq  is  the  intensity  of  the  light  incident  on  the  sample  and  /  is  the  intensity  of  the  light 
that  returns  from  the  sample  and  is  measured  by  the  instrument.  The  sample  absorbance 
(A)  is 

A  =  -log,„(fl)  (2) 

and  is  proportional  to  analyte  concentration  (c)  under  certain  conditions.  The  Beer- 
Lambert  law  is 

A  =  abc  (4) 

where  b  is  the  path  length,  the  distance  the  radiation  traveled  through  the  sample,  and 
a  is  the  molar  absorptivity,  a  constant  characteristic  of  the  analyte.  For  a  more  complex 
sample,  the  sample  absorbance  (Ajj)  is  a  linear  sum  of  the  absorbances  of  all 
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components. 

Aj  ~  k^Ajk  “  (k^^jk  ^ik)  (1) 

It  is  assumed  that  all  of  the  radiation  travels  the  same  distance  through  the  sample  so  that 
b;  is  constant  for  all  rays.  For  diffuse  reflectance,  this  is  not  strictly  true,  but  it  is  often  a 
useful  approximation. 

Other  common  transformations  include  mean  centering,  baseline  subtraction,  and 
derivatives.  The  equation  for  mean  centering  is 

)mean  centered  ~  ^\  *  ^(^ij)  (5) 

where  is  the  column  vector  of  log(1  /R)  for  the  wavelength  and  all  samples.  Baseline 
subtraction  is 

(^i)baseline  corrected  ~  *  ^i,Jb  (®) 

where  r,  is  the  spectrum  for  sample  /  and  is  the  absorbance  at  the  Jt,*”  wavelength. 


An  estimation  of  the  second  derivative  using  a  finite  difference  function  is 

f’  *  ^  +  y  /  9^  (7) 

where  fj  is  the  intensity  at  the  f'  wavelength  and  g  represents  the  gap  size  of  the 
derivative. 

Calibration 

Calibration  is  used  for  continuous  constituent  values,  such  as  chemical 
concentrations  or  physical  properties.  Calibration  models  must  be  validated  before  they 
may  be  used  for  prediction.  In  calibration,  a  model  is  developed  that  relates  the  spectra 
data  matrix  (R),  composed  of  the  spectra  vectors  fj,  to  the  constituent  matrix  (C), 
composed  of  the  elements  Cj,,,  for  a  set  of  samples,  called  the  training  set  or  the 
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calibration  set  (designated  by  the  subscript  c)-  A  second  set  is  used  for  validation.  The 
model  is  used  to  predict  the  constituent  values  of  the  second  set.  Poor  prediction 
indicates  that  either  irrelevant  data,  such  as  instrument  noise  was  included  in  the 
calibration  model  (overfitting)  or  that  there  is  variance  in  the  validation  set  that  was  not 
present  in  the  training  set,  such  as  a  new  component.  Once  the  model  has  been 
validated,  it  may  be  used  to  predict  the  constituents  of  unknown  samples  (designated  by 
the  subscript  p). 

Another  method  of  validation,  leave-one-out  cross  validation,  is  useful  when  there 
are  only  a  few  spectra  available.  Instead  of  dividing  the  set  in  two,  half  for  the  calibration 
set  and  half  for  the  prediction  set,  all  samples  except  one  are  used  for  the  calibration  set 
and  the  remaining  sample  becomes  the  prediction  set.  This  is  repeated  /  times,  leaving 
out  and  predicting  each  spectra  once.  The  result  gives  an  estimate  of  the  accuracy  of 
the  model  that  uses  all  /  samples.  Since  the  spectra  acquired  from  different  sites  on  the 
same  wound  may  be  highly  correlated,  this  approach  was  modified:  all  spectra  from  a 
single  patient  were  left  out  and  predicted  together. 

Classical  Least  Squares 

The  calibration  model  generally  consists  of  a  vector  of  matrix  of  coefficients  (S)  that 
relates  the  constituent  matrix  to  the  spectral  matrix.  The  coefficients  are  typically 
estimated  by  least  squares  or  similar  methods.  The  most  straight  forward  methods 
calculate  S  assuming  the  Beer-Lambert  relationship. 

Classical  least  squares  (CLS)  is  a  direct  calibration  technique  based  on  the  Beer- 
Lambert  law,  which  in  matrix  form  is 

R  =  C  S  +  Ef,  (8) 

where  S  is  the  estimated  pure  component  spectra  (Sj^  «  aj^  bjj)  and  Er  is  the  matrix  of 
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spectral  errors.  The  S  matrix  is  determined  by 

S  =  Re  Ce"  (9) 

where  Rg  and  Cg  are  the  spectral  and  concentration  matrices  of  the  training  set.  is 
the  Moor-Penrose  pseudolnverse:^ 

C*  =  (C^  C)-’  (10) 

By  solving  these  equations,  Er  is  minimized  in  a  least  squares  sense.  The  estimated 
concentrations  (Cp^)j  of  all  k  components  in  the  unknown  sample  with  spectrum  (fp^); 
is 

=  (rp')i  S'  (11) 

It  is  often  useful  to  include  constraints,  such  as  requiring  all  elements  in  S  and  (Cp^)i  to 

be  greater  than  or  equal  to  zero. 

Multiple  Linear  Regression 

When  the  sample  spectra  cannot  be  approximated  as  a  simple  sum  of  pure 
component  spectra,  either  because  some  of  the  components  are  unknown  or  because 
Beer’s  law  is  inapplicable,  indirect  calibration  in  used.  The  basic  model  is 

C  =  R  P  +  Ec  (12) 

where  Eq  is  the  error  in  the  concentration  matrix  and  P  is  the  /  by  J  matrix  of  coefficients, 
similar  to  S.  C  may  be  any  constituents,  not  just  chemical  concentrations.  The  equation 
is  solved  to  minimize  Eq,  and  prediction  is 

Since  estimating  P  involves  calculating  the  pseudoinverse  of  Rg,  J  must  be  small.  If  the 
spectra  originally  contain  a  large  number  of  wavelengths,  there  are  two  methods  of 
simplifying  the  spectral  data:  selecting  a  few  relevant  wavelengths  and  factoring  the 
whole  spectrum. 


17 


Multiple  linear  regression  (MLR),  also  known  as  inverse  least  squares^  and  the 
P*matrix  approach^’^,  is  an  inverse  least  squares  model  that  uses  the  spectral  data  from 
a  small,  selected  set  of  wavelengths.  The  model  is 

=  R  p/  +  (14) 

where  pj  is  the  vector  of  regression  coefficients  for  the  constituent  and  is  the  error 
in  the  concentration  matrix. 

The  most  common  method  for  choosing  the  wavelengths  to  use  in  the  MLR  model 
is  to  calculate  the  correlation  coefficients  between  all  wavelengths  and  the  constituent  of 
interest.  The  wavelength  with  the  greatest  correlation  in  chosen  as  the  first  wavelength 
for  the  model.  The  variance  associated  with  the  correlation  is  removed  from  the  spectral 
and  constituent  data  and  a  second  wavelength  is  chosen,  and  so  on.  Other  methods 
also  exist. 

Principal  Component  Regression 

Principal  component  regression  (PCR)^’^  uses  the  same  fundamental  definition 
as  MLR  but  decomposes  the  full  spectral  matrix  R  into  two  factor  matrices  T  and  P, 

R  =  T  P  +  Er  (15) 

which  can  be  calculated  by  singular  value  decomposition: 

R  =  V  S  +  Ep  (16) 

where  V  is  /  by  >A,  S  \sA  by  A,  and  is  A  by  J,  and  A  is  the  number  of  factors  (principal 
components).  T  and  P  are  VS  and  U^,  respectively.  The  matrix  T  contains  "scores," 
linear  combinations  of  concentratioadata,  and  P  contains  "loadings,"  linear  combinations 
of  spectral  data.  Scores  and  loadings  are  difficult  to  interpret  in  terms  of  the  individual 
components.  The  concentrations  can  be  estimated  from  the  scores. 


c  =  T  q  + 


(17) 
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q  »  t/  c.  (18) 

Cp^Rp^Vc,  (19) 

where  q  is  the  vector  of  coefficients  that  relates  c  to  T.  If  Cp  is  not  calculated,  this  method 
is  known  as  principal  component  analysis  (PCA)  and  does  not  require  any  constituent 
values.  This  method  is  useful  as  a  classification  technique,  particularly  when  constituent 
values  are  unavailable. 

Classification 

When  constituent  data  is  discrete,  such  as  the  place  of  origin  of  the  samples  or 
"acceptable"  and  "unacceptable",  classification  algorithms  are  used.  The  most  common 
is  K  Nearest  Neighbors  (KNN).  The  samples  are  treated  as  points  in  a  multidimensional 
space  and  their  "nearness"  to  each  other  are  calculated: 

0,j  =  lir,  -  fjil  (20) 

where  D,j  is  the  distance  between  samples  i  and  /,  and  fj  and  fj  are  the  spectral  or  score 
vectors  for  the  two  samples.  It  is  assumed  that  samples  within  a  class  will  be  closer  to 
each  other  than  to  samples  in  another  class. 

Like  calibration,  KNN  classification  requires  a  training  set  of  samples  whose 
classifications  are  known.  However,  there  is  no  matrix  of  regression  coefficients,  an 
unknown  sample  is  simply  compared  to  all  samples  in  the  training  set,  one  at  a  time.  If 
K  =  1 ,  the  unknown  sample  is  assigned  to  the  class  of  the  training  set  sample  that  it  is 
closest  to.  Otherwise  the  classifications  of  the  K  nearest  training  samples  are  considered. 

The  data  used  for  KNN  is  often  scaled  (weighted)  so  that  the  variables  which  have 
the  largest  correlation  to  the  classification  are  bigger  than  those  which  are  less  relevant. 
One  common  method  of  determining  which  variables  are  relevant  is  the  variance 
weight,^  which  is  the  ratio  of  the  intercategory  variances  to  the  sum  of  the  intracategory 
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variances.  The  equation  is 

=  2  *  KVWJa."  +  (VWb)ab“  -  (2/N.NJS(.stbl  / 

[(1/N.)sK-x.,J^  +  (1/A)jE0t,-x^,J^l  (21) 

where  and  A/^,  are  the  number  of  samples  in  classes  a  and  b,  x  are  the  values  of  the 
variable  in  question,  and  x,  ^  is  the  mean  of  x  for  samples  in  class  a. 

RESULTS 
Case  Study 

A  case  study  is  included  here  to  illustrate  typical  results  for  shallow  and  full 
thickness  burns.  Patient  #14  had  a  small  scald  burn  on  her  foot.  Hot  water  had  run 
down  into  her  high-top  boot,  producing  a  shallow  partial  thickness  burn  at  the  top  and 
a  deep  burn  where  the  water  pooled  at  the  bottom.  In  the  spectrum  of  the  shallow  partial 
thickness  bum  (Rgure  15a)  the  oxy-  and  deoxy-  hemoglobin  bands  (400-600  nm)  are  very 
intense,  because  the  epidermis  was  absent  and  the  wound  was  inflamed.  In  the  near 
infrared  region  (700-1800  nm)  there  are  three  absorbance  bands  that  arise  from  the 
overtones  of  0-H  stretching  modes  of  water  (960, 1150,  and  1450  nm).  The  discontinuity 
in  the  spectrum  at  900  nm  is  an  instrumental  artifact,  caused  by  a  filter  change. 

The  spectrum  of  the  deep  burn  (Figure  15b)  was  markedly  different.  The  baseline 
was  considerably  higher  and  the  water  bands  in  the  near  infrared  were  less  intense.  The 
cause  of  these  differences  is  unknown.  The  560  nm  hemoglobin  bands  are  also  less 
intense,  as  there  is  less  blood  in  the  injured  skin.  All  of  the  above  had  been  expected. 
At  630  nm,  however,  the  is  an  unanticipated  band.  The  intensity  of  this  band,  evident  as 
a  shoulder  on  the  oxy-/deoxy-  hemoglobin  peak,  suggested  that  it  was  due  to  another 
form  of  hemoglobin.  The  band  was  attributed  to  the  acid  form  of  methemoglobin  which 
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has  a  band  at  630  nm  (Figure  1a).  This  compound  is  relatively  well  known  and  is  present 
in  normal  subjects  at  the  0.5-1. 0%  level.  Other  hemoglobin  derivatives  may  also  be 
present  in  smaller  quantities. 

Principal  Component  Models 
Raw  Spectra 

The  spectra  from  all  of  the  burn  wounds  are  shown  in  Figure  16.  Twelve  spectra, 
acquired  from  discolored  (black,  brown,  or  green;  n  =  10)  and  first  degree  (n=2)  burns 
contain  spectral  features  that  are  not  relevant  to  the  discrimination  of  shallow  and  deep 
burns.  They  were  not  included  in  any  of  the  calibration  models,  to  avoid  making  the 
models  more  complicated  than  necessary,  but  were  included  in  predictions  to 
demonstrate  that  they  would  not  be  misclassified.  These  spectra  are  shown  in  Figure  16 
and  subsequent  figures  as  dotted  lines. 

The  spectral  regions  which  correlate  to  burn  depth  were  located  by  calculating  the 
variance  weights  (Equation  21)  for  the  absorbances  at  each  wavelength  (Figure  17). 
There  is  some  correlation  between  550  and  580  nm,  but  little  or  none  (variance  weight 
=  2)  anywhere  else.  This  indicates  that  the  hemoglobin  absorbances  (550-580  nm)  have 
some  correlation  to  burn  depth,  which  is  in  accord  with  one  of  the  fundamental  tenets  of 
clinical  assessment,  that  the  darker  pink  a  wound  is,  the  more  likely  it  is  to  be  shallow. 
The  low  correlation  in  other  portions  of  the  spectra  indicates  that  the  baseline  offsets  are 
unrelated  to  burn  depth  and  are  probably  caused  by  irreproducibility  in  the  source  / 
sample  /  detector  geometry. 

Since  the  major  spectral  features  do  not  correlate  strongly  to  burn  depth,  the 
spectra  were  analyzed  by  Principal  Component  Analysis  (Equations  15-19)  to  locate  minor 
features  which  might  be  more  relevant.  The  spectra  would  be  expected  to  display  at  least 
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four  major  types  of  variation,  corresponding  to  baseline  offsets,  multiplicative  factors  (path 
length  and  total  blood  concentration),  oxygenation  of  the  blood,  and  methemoglobin 
concentration.  Minor  variations,  representing  minor  chemical  components  or  spectral 
non-linearities  would  also  be  expected.  The  first  four  principal  component  loadings 
(spectral  domain.  Figure  18a,b)  show  that  there  are  indeed  at  least  four  significant  types 
of  orthogonal  spectral  variance.  The  next  four  (Figure  18c,d)  model  only  a  small  amount 
of  the  total  variance.  It  is  possible,  but  unlikely  that  these  later  factors  are  significant. 

The  principal  component  scores  (concentration  domain.  Figure  19)  show  that  the 
first  two  factors,  are  dominated  by  irrelevant  effects  of  sample/instrument  geometry. 
Similar  effects  have  been  observed  in  the  diffuse  reflectance  spectra  of  other  irregularly- 
shaped  samples.^’  The  third  and  fourth  principal  components,  however,  are  correlated 
to  burn  depth.  The  variance  weights  for  the  first  eight  principal  components  are  2.3,  2.5, 
2.7,  3.2,  2.5,  2.3,  2.0,  and  2.0. 

The  distinction  between  "shallow"  and  "deep  is  a  somewhat  arbitrary  division  of  a 
continuous  variable,  the  fraction  of  dermis  destroyed.  It  might  be  expected  that  in  a 
multidimensional  space,  the  scores  for  shallow  burns  might  be  located  on  one  side  and 
those  for  deep  burns  on  the  other,  with  the  borderline  injuries  inbetween  (Figure  19b). 
Shallow  burns  would  most  likely  be  closest  to  other  shallow  burns,  etc.,  so  KNN  should 
give  reasonable  results.  The  KNN  algorithm  in  Pirouette^^  successfully  classified  1(X) 
(91%)  of  the  110  spectra,  using  4  principal  components.  The  optimal  number  of 
neighbors  to  consider  was  found  to  be  2.  Cross  validation,  leaving  out  one  patient’s 
spectra  each  time  (MATLAB  program  written  by  the  author),  successfully  classified  101 
(92%),  using  6  principal  components  and  considering  3  neighbors. 

Assuming  that  there  is  only  one  continuous  variable  which  relates  to  burn  depth, 
another  possible  method  of  classification  is  to  define  a  multidimensional  plane  between 
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the  two  classes.  All  scores  on  one  side  of  the  plane  should  represent  shallow  burns,  and 
all  scores  on  the  other  side  should  represent  deep  injuries.  This  plane  may  be  found  by 
rotating  the  axes  so  that  the  first  axis  corresponds  to  the  direction  of  the  variance  in  the 
data  which  most  clearly  distinguishes  between  the  two  classes.  The  discrimination  plane 
is  then  perpendicular  to  the  first  rotated  axis. 

The  algorithm  used  for  single-plane  discrimination  involved  5  steps.  The  first  step 
was  to  mean-center^  the  scores.  Then  the  mean-centered  scores  were  range  scaled^ 
and  feature  weighted,^  according  to  the  variance  weight.  The  scaled  data  was  then 
rotated,  two  dimensions  at  a  time,  starting  with  those  which  have  the  greatest  variance 
weights.  The  best  rotation  was  defined  as  that  which  produced  the  largest  variance 
weight  for  the  first  of  the  two  dimensions.  Once  the  best  rotation  was  found,  the 
discrimination  plane  was  defined  perpendicular  to  the  first  axis,  located  so  that  the 
number  of  misclassified  calibration  samples  was  minimized.  The  samples  were  classified 
according  to  which  side  of  the  discrimination  plane  their  rotated  scores  were  located.  The 
MATLAB  code  for  these  programs  are  in  Appendix  B. 

The  scores,  rotated  in  the  first  four  dimensions,  are  shown  in  Figure  20.  The 
discrimination  plane  is  also  shown,  projected  as  a  line.  All  but  7  of  the  samples  are 
correctly  classified  by  this  plane.  Including  more  than  four  dimensions  does  not  improve 
classification,  since  only  the  first  four  principal  components  contain  much  relevant 
information.  Cross  validation,  using  4  principal  components  correctly  classified  102 
(93%).  The  single-plane  discrimination  Is  slightly  better  than  KNN  because  only  one 
dimension  contains  relevant  information.  The  model  is  simpler  and  less  influenced  by  any 
single  calibration  sample. 

The  loadings,  rotated  in  the  same  manner  (Figure  21),  show  that  there  is  a  strong 
positive  correlation  between  the  absorbance  at  630  nm  and  the  depth  of  the  burn.  This 
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is  probably  due  to  the  absorbance  of  methemoglobin  which  appears  to  be  characteristic 
of  deep  burns.  There  is  also  a  strong  negative  correlation  at  570  nm,  where 
oxyhemoglobin  absorbs.  It  is  reasonable  to  assume  that  the  more  healthy,  oxygenated 
hemoglobin  a  wound  contains,  the  more  likely  it  is  to  heal  quickly.  The  first  rotated 
loading  is  small  for  the  wavelengths  650-700  nm  indicating  that  this  region,  which  contains 
little  hemoglobin  absorbance,  is  not  useful  for  prediction.  This  further  supports  the 
previous  results  indicating  that  the  baseline  offsets  are  unrelated  to  burn  depth. 
Subsequent  rotated  loadings  (not  shown)  are  not  correlated  to  burn  depth  (variance 
weights  <  2.2). 
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Baseline  subtracted  spectra 

Since  the  baseline  offsets  are  useless,  they  can  be  removed  without  losing 
information,  approximating  the  baseline  by  the  apparent  absorbance  at  700  nm  (Equation 
6,  Jb  -  700  nm).  Among  the  resulting  spectra  (Figure  22)  the  spectral  differences 
between  the  two  classes  are  more  obvious.  The  intensity  of  the  hemoglobin  absorbance 
(520-580  nm)  is,  on  the  average,  greater  for  shallow  injuries  than  for  deep,  but  this 
variable  alone  is  insufficient  for  diagnosis.  Considering  only  the  intensity  of  the 
absorbance,  what  the  eye  sees,  only  69%  of  the  injuries,  at  best,  could  be  correctly 
classified.  Clinical  assessment  considers  a  few  other  variables  and  is  potentially  more 
accurate. 

Since  the  variance  due  to  the  offsets  has  been  removed,  fewer  principal 
components  should  be  needed  to  describe  the  data  set.  The  first  four  loadings  (Figure 
23)  are  similar  to  those  from  the  raw  spectra  (Figure  18).  The  subsequent  loadings, 
however,  contain  more  noise.  The  second  and  third  principal  components  are  most 
highly  correlated  (Figure  24),  rather  than  the  third  and  fourth  (Figure  19),  and  the  variance 
weights  for  the  principal  components  are  greater,  2.5,  4.6,  2.9,  and  2.2  for  the  first  four. 
Prediction  results  from  KNN  and  single-plane  discrimination  are  not  significantly  better, 
but  the  model  is  simpler. 

Normalized  spectra 

Although  the  amount  of  blood  visible  in  a  wound  is  a  good  indication  of  its  viability, 
it  is  an  unreliable  indicator.  Edema  (inflammation)  and  hematoma  (bruising)  both  increase 
the  amount  of  blood  without  increasing  the  healing  potential.  In  this  study,  the  diagnostic 
errors  made  by  the  physicians  were  caused  either  by  injuries  that  were  outside  the  range 
of  their  experience,  such  as  an  unusual  chemical  burn,  or  by  deep  injuries  that  appeared 
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pink  instead  of  white.  To  correctly  diagnose  these  wounds  it  is  helpful  to  focus  on  the 
second  factor  which  relates  to  burn  depth,  the  methemoglobin  content. 

To  remove  the  variance  caused  by  variable  amounts  of  blood  in  the  skin,  the 
spectra  were  normalized  to  the  maximum  of  the  average  spectrum  (570  nm  was  chosen). 
The  normalization  equation  is 

(^i) normalized  ~  (^i)baseline  subtracted  /  ^i,570nm  (22) 

where  r^sronm  'S  the  apparent  absorbance  of  the  /*”  spectrum  at  570  nm.  The 
methemoglobin  absorbance  at  630  nm  is  clearly  visible  in  the  normalized  spectra  (Figure 
25). 

The  removal  of  the  multiplicative  factors  is  inaccurate  since  this  normalization  is 
somewhat  simplistic,  so  the  rank  of  the  spectral  data  set  does  not  decrease.  But  since 
the  data  is  simplified,  the  correlation  to  burn  depth  increases  (Figure  26).  The  loadings 
(Figure  27)  are  similar  to  those  of  the  baseline-corrected  spectra  (Figure  23),  but  the 
scores  of  the  first  two  principal  components  (Figure  28)  are  more  highly  correlated.  The 
variance  weights  for  the  first  four  principal  components  are  3.5,  5.3,  2.8,  and  2.2. 
Prediction  results  are  nearly  identical  to  those  from  the  baseline-corrected  data. 

Hemoglobin  Spectra  Model 
Single  Layer  Model 

The  above  results  suggest  that  there  is  a  single  spectral  feature  or  single 
combination  of  interrelated  features  which  is  indicative  of  the  depth  of  a  burn  wound. 
Further,  in  the  normalized  spectra,  the  630  nm  absorbance,  attributed  to  the  acid  form  of 
methemoglobin,  correlated  to  burn  depth.  This  suggests  that  the  concentration  of  acidic 
methemoglobin  in  the  skin,  relative  to  the  total  blood  volume,  is  proportional  to  the 
amount  of  tissue  damage  and,  hence,  to  the  length  of  time  required  for  healing. 
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Unfortunately,  there  is  no  way  to  test  this  hypothesis  directly.  Even  if  the  patients 
consented  to  having  skin  samples  removed  for  analysis,  conventional  methods  of 
hemoglobin  speciation  are  not  applicable  to  solid  tissue  samples. 

One  alternative  is  to  attempt  to  reproduce  the  observed  reflectance  spectra, 
starting  from  the  spectra  of  pure  hemoglobins.  Assuming  that  (a)  the  particles  in  the  skin 
scatter  light  isotropically,  (b)  the  hemoglobin  is  uniformly  distributed,  (c)  there  is  no 
specular  reflectance  or  stray  light,  and  (d)  the  only  absorbing  compounds  in  the  skin  are 
hemoglobins,  then  the  reflectance  spectra  of  the  skin  could  be  modeled  as  a  linear  sum 
of  pure  hemoglobin  spectra,  using  Classical  Least  Squares  (Equations  8-11).  None  of  the 
assumptions  are  actually  true  and  the  hemoglobins  used  to  acquire  the  basis  spectra 
were  not  100%  pure,  so  the  model  is  not  expected  to  fit  the  data  precisely.  The  spectra 
of  the  hemoglobins  (the  matrix  S)  are  shown  in  Figure  29,  including  oxyhemoglobin, 
deoxyhemoglobin,  acid  methemoglobin,  base  methemoglobin,  and  a  baseline.  The 
spectra  were  fit  to  a  reflectance  spectrum  of  a  burn  wound,  and  a  new  spectrum  was 
reconstructed  from  the  estimated  hemoglobin  concentrations.  This  was  repeated  for  all 
burn  spectra. 

Two  examples  are  shown  in  Figure  30.  The  reconstructed  spectra  are  not  exactly 
the  same  as  the  measured  spectra,  since  the  assumptions  are  incorrect.  This  simplistic 
model  does,  however,  give  an  estimate  of  methemoglobin  in  a  burn,  information  which 
has  previously  been  unavailable.  The  resulting  predictions  are  less  accurate  (79%)  than 
those  from  the  empirical  PCA  models,  indicating  that  the  methemoglobin  concentration 
alone  is  insufficient  to  predict  healing  or  that  the  fitting  errors  were  too  large,  probably 
both.  However,  the  success  of  this  model  in  discrimination  does  strongly  support  the 
hypothesis  that  oximetric  information  correlating  to  burn  depth  can  be  obtained  from  the 
spectra. 


27 


Two  Layer  Model 

Of  the  assumptions  used  in  the  above  model,  two  are  approximately  true,  one 
might  be  true,  and  one  is  definitely  false.  The  hemoglobin  is  not  distributed 
homogeneously;  in  most  cases  there  is  a  layer  of  bloodless,  dead  tissue  (eschar)  at  the 
surface.  This  layer  reflects  incident  light,  reducing  the  amount  that  interacts  with  the 
blood  beneath.  The  reflectance  from  a  non-absorbing  layer  is  dependent  on  light 
scattering  and  is  a  function  of  wavelength. 

If  the  amount  of  light  of  wavelength  /,  reflected  from  the  sample  Is  /j,  and  the  total 
incident  light  is  /qj,  then  the  instrument  response  (/?obs,j)  is 

«obs.i  =  h  /  /o.i  (23) 

«obs.i  =  iUi  +  g  /  (/soj  +  g  (24) 

where  /,j  is  the  intensity  reflected  from  the  non-absorbing  eschar  layer,  /so,j  is  the  intensity 
of  the  light  that  passes  through  the  first  layer  and  enters  the  second,  blood-filled  layer, 
and  is  the  intensity  of  the  light  returning  from  the  second  layer.  The  fraction  of  the 
incident  radiation  that  would  be  diffusely  reflected  from  the  second  layer,  if  there  were  no 


specular  reflectance  and  if  the  first  layer  were  removed,  is 

«».i  =  'si  /  Ui  (25) 

and  is  related  to  the  absorbance  of  the  blood  in  the  second  layer  (Aj  j), 

^.i«iog(1/«s,i)  (26) 

It  can  be  estimated  from  the  instrument  response; 

«s.i  =  (/oj  «0bs.i  -  /r.i)  /  (/o.i  -  /r.i)  (27) 

«..i  =  (7obs.i  -  S)  /  (1  -  Sj)  (28) 

Sj  =  l,i  /  /o,i  (29) 

if  Sj  is  known. 


For  the  burn  spectra,  S  is  not  known.  Instead,  S  was  approximated  by  a  second 
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order  polynomial,  a  function  of  the  wavelength  (A): 

S»s„/(A»)  +  +  S:f(X‘)  (30) 

/?3  should  be,  approximately,  the  product  of  the  transmittances  (TJ  of  the  hemoglobin 
species  {k)  times  their  concentrations  (cj. 

«s  «  wH  K  T,)  (31) 

And  therefore 

(7-obs.i  -  S)  /  (1  -  Sj)  «  ,n  (c,  W  (32) 

There  are  3+/f  variables  which  must  be  solved  simultaneously  (Sq,  s,,  S2,  c,,  Cg, ...  cj  by 
using  a  non-linear  curve  fitting  algorithm  and  minimizing  the  spectral  fitting  error: 

error  =  (To^sj  -  Sj)  /  (1  -  Sj)  -  ^k)  (33) 

The  Gauss-Newton  algorithm  supplied  by  Mathworks,^^  gave  satisfactory  results,  under 
the  following  constraints: 

Ck  ^  0  (34) 

0  ^  S,  ^  1  (35) 

The  reconstructed  spectrum  (A^bs)  is  then 

7-ob,.i  =  7-3.j  (l-Sj)  +  Sj  (36) 

>^obs.j « -  iog(rob,.i)  (37) 


The  transmission  spectra  of  hemoglobins  and  the  functions  /(A°),  /(A’),  and  /(A^) 
are  shown  in  Figure  31.  The  reconstructed  spectra  for  the  two  examples  are  shown  in 
Figures  32  and  33.  In  general,  the  fit  was  not  improved  if  S  was  assumed  to  be  0**’  order 
with  respect  to  wavelength  (s,  =  0  and  S2  =  0).  When  S  was  1®*  order,  there  was 
significant  improvement;  nearly  all  of  the  systematic  distortion  was  accounted  for.  Adding 
the  2"*^  order  improved  some  of  the  more  unusual  spectra,  but  generally  had  little  effect 
on  the  "typical"  spectra.  The  1®*  order  model  indicated  that  deep  burns  usually  have  more 
than  9%  acid  methemoglobin  (Figure  34)  and  correctly  classified  90%  of  the  injuries. 
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Skin  Temperature  Model 

In  the  second  year,  spectra  were  acquired  from  the  forearms  of  healthy  volunteers. 
Their  skin  temperature  was  varied  from  9  to  35  °C  and  measured  with  a  thermistor  at  the 
surface.  The  spectra,  after  second  derivative  transformation,  were  analyzed  by  Multiple 
Linear  Regression  (Equations  12-14)  and  were  found  to  correlate  to  the  skin  temperature 
(correlation  coefficient  R  =  0.982),  with  a  standard  error  of  prediction  of  1.4  °C.  This 
experiment  is  described  in  the  previous  report.'*®  The  equation  was 

(  =  22.7  +  (1.21X10«)  r-,„^  -  (1.44x10“)  ■  (3.85x10=)  (38) 

where  r"ioi5nm  'S  the  second  derivative  (Equation  7,  gap =23)  at  1015  nm  and  t  is  the 
temperature  in  ®C. 

The  model  for  estimating  skin  temperature  in  normal  skin  was  used  to  estimate  the 
temperature  of  the  burn  wounds  (Figure  35).  The  results  suggest  that  deeper  burns  are 
cooler  than  shallow  burns  (Figure  36).  The  average  estimated  temperatures  were  32  ± 
7  and  27  ±  7  for  shallow  and  deep  burns  (P<0.01).  Thermography,  using  emitted 
infrared  light,  has  also  shown  that  deep  burns  are  cooler  than  shallow,,,  ,2  but  that  such 
measurements  are  extremely  sensitive  to  cooling  caused  by  evaporation  of  surface  water. 
NIR  reflectance  should  be  less  sensitive  to  surface  effects,  since  it  measures  that 
temperature  within  the  skin  and  below  it.  The  model  used  here  is  not  optimal;  it  is  derived 
from  healthy  skin.  A  more  accurate  model  could  be  obtained  empirically,  if  the 
temperatures  of  wounds  could  be  measured  by  a  contacting  or  invasive  probe,  or 
theoretically,  if  the  spectral  non-linearities  were  better  understood. 

Although  several  features  in  the  NIR  spectra,  including  water  temperature  correlated  to 
burn  depth,  none  were  as  useful  as  the  hemoglobin  bands,  nor  did  they  improve 
prediction.  It  is  likely  that  there  is  still  much  useful  information  in  these  spectra,  but  in  the 
absence  of  physiological  reference  data,  this  information  could  only  be  extracted  by 


30 


rigorous  application  of  diffuse  reflectance  theory. 

BDI  and  IBDI  Simulations 

Now  it  is  possible  to  explain  the  predictive  ability  of  the  old  Burn  Depth  Indicators 
(imaging  and  non-imaging).  Both  Instruments  used  baseline-subtracted,  single¬ 
wavelength  normalized  data.  The  baseline  offsets  were  approximated  by  the  instrument 
response  of  the  NIR  LED  (880  nm)  or  filter  (880-1100  nm)  and  the  remaining  data  was 
normalized  to  the  response  of  the  Green  (560  nm).  The  methemoglobin  concentration 
was  then,  inadvertently,  approximated  by  the  normalized  response  of  the  Red  (640  nm). 
The  exact  mathematics  are  slightly  different,  because  instruments  used  the  instrument 
response  (R)  rather  that  log(1/R),  but  the  principle  is  the  same. 

In  order  to  compare  the  prediction  results  from  our  spectrophotometer  with 
predictions  from  the  Burn  Depth  Indicator  (BDI)  and  the  Imaging  Burn  Depth  Indicator 
(IBDI)  built  previously,  the  BDI  and  IBDI  responses  were  simulated  from  our  spectra.  The 
emission  spectra  of  the  light  emitting  diodes  (LED’s)  in  the  BDI  were  simulated  as 
Gaussian-shaped  peaks  with  maxima  at  550,  640,  and  880  nm  and  band  widths  (fwhm) 
of  70  nm.^’’“  The  spectra  of  the  IBDI  filters  (Figure  37b)  were  obtained  from  the 
literature.^  The  responses  for  the  near  infrared  LED  and  filter  are  somewhat  uncertain, 
because  our  spectral  data  between  750  and  890  nrn  is  unreliable  and  was  not  used.  The 
simulated  results  (Figures  38a  and  38b),  however,  are  reasonably  close  to  those  obtained 
in  previous  work.’®'“  For  the  few  wounds  for  which  IBDI  images  were  acquired,  the 
simulated  results  agree  with  the  images.  The  burn  sites  used  in  this  study  seem  to  have 
a  slightly  wider  range,  both  for  the  Red/NIR  and  Green/NIR  ratios,  than  those  used  in 
previous  studies,  which  is  consistent  with  our  attempt  to  select  the  widest  variety  of 
injuries.  The  accuracy  of  the  BDI  and  IBDI  simulations  for  prediction  burn  depth  was 
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about  77  and  79%,  respectively.  VanLiew  and  Moore  reported  79  and  84%  accuracies 
for  different  sets  of  indeterminate  burns  for  the  two  instruments.^®'^ 

Proposed  Model  for  New  Imaging-Bum-Depth-Indicator 

These  results  show  that  the  model  used  for  the  Burn  Depth  Indicators  was 
essentially  sound  and  from  the  perspective  of  burn  physiology,  basically  correct,  though 
its  explanation  was  erroneous.  The  choice  of  wavelengths,  however,  could  be 
substantially  improved.  Rrst,  the  filter  originally  selected  for  baseline  subtraction  (880- 
1100  nm)  measures  a  large  absorbance  band  at  970  nm  (Figure  37)  as  well  as  the 
offsets.  This  band  arises,  in  part,  from  the  water  in  and  on  the  skin.  The  water 
absorption  is  variable  but  not  well  correlated  to  burn  depth.  The  880-1100  nm  filter  also 
measures  the  absorbances  of  hemoglobins  (Figure.  1b),  but  the  path  length  in  this  region 
is  much  longer  than  in  the  visible  region.  It  is  uncertain  whether  the  blood  in  the  volume 
sampled  is  representative  of  the  blood  in  the  skin  or  primarily  below  the  skin. 
Furthermore,  the  BDI  and  IBDI  simulations  produced  comparable  results,  in  spite  of  the 
fact  that  the  wavelength  ranges  used  for  baseline  correction  were  significantly  different. 
Therefore,  little  relevant  information  is  lost  and  some  interference  is  avoided  by  moving 
the  baseline-correction  filter  to  650  or  700  nm,  where  the  absorbances  are  significantly 
lower. 

Second,  the  normalization  filter  could  be  narrowed  and  shifted  to  produce  a  more 
reliable  response.  There  is  a  pronounced,  non-linear  distortion  of  the  spectra  which 
affects  the  550-580  nm  region  where  the  hemoglobin  absorbance  is  particularly  strong. 
In  the  spectra  of  whole  blood  (Figure  la),  the  two  oxyhemoglobin  absorption  bands  (535 
and  570  nm)  have  similar  intensities.  In  the  burn  spectra,  however,  the  535  nm  band  is 
decreased  relative  to  the  570  nm  band,  and  the  degree  of  distortion  is  unrelated  to  burn 


depth.  The  distortion  is  probably  caused  by  scattering,  and  scattering  is  highly 
wavelength  dependent.  Therefore,  the  scattering  at  630  nm  is  likely  to  be  more  similar 
to  that  at  570  nm  than  535  nm,  and  the  former  should  give  a  more  accurate  normalization 
for  the  630  nm  methemoglobin  band.  Simulations  of  fixed-filter  instruments  using  various 
possible  normalization  filters  confirm  this.  The  normalization  LED  for  the  BDI  had  a  80  nm 
band  width  and  was  at  560  nm,  sampling  both  oxyhemoglobin  bands  equally.  The  filter 
in  the  IBDI  at  540  nm  was  worse,  sampling  only  the  most  distorted  portion  of  the 
oxyhemoglobin  spectrum.  A  filter  with  a  narrower  band  width,  30  nm  or  less,  and  located 
at  570  nm  should  give  a  more  reliable  normalization. 

Third,  the  methemoglobin  filter  should  be  optimized.  The  filter  used  previously  had 
a  80  nm  band  width  and  was  centered  at  640  nm.  The  methemoglobin  absorbance  peak 
is  at  630  nm,  and  for  the  normalized  spectra,  the  wavelength  that  correlates  most  to  burn 
depth  is  625  nm  (Figure  26).  The  band  width  of  the  methemoglobin  band  is 
approximately  40  nm.  A  narrower  filter,  40  nm  fwhm  or  less,  centered  at  630  nm  should 
give  more  precise  estimates  of  methemoglobin  absorbance. 

One  possible  filter  selection  might  include  three  20-nm  band  pass  filters  at  650, 
570,  and  630  nm.  If  the  apparent  absorbances  (log  (1/R))  of  the  three  filters  are  Agg^, 
Ajtq,  and  Aggo  respectively  (Figure  39),  then  the  burn  depth  (D)  is  estimated  by 

^  “  (^630  ■  ^65o)  /  (^570  '  ^65o)  (^) 

The  discrimination  point,  calculated  the  same  way  as  the  discrimination  plane  for  the 
rotated  PCA  scores,  was  found  to  be  0.104.  Injuries  for  which  D  is  greater  than  0.104 
were  deep  and  the  rest  were  shallow.  One  hundred  one  (92%)  of  the  sites  were  correctly 
classified  by  this  model,  and  97  (88%)  were  classified  during  cross  validation.  This  model 
was  significantly  more  accurate  (88%)  than  the  IBDI  simulation  (79%)  for  the  104  sites 
simulated  (P<0.06  for  paired  data).  Small  variations  in  the  band  width  or  mean 
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wavelength  of  the  filters  do  not  significantly  affect  prediction  ability,  indicating  that  an 
instrument  using  this  algorithm  would  be  relatively  tolerant  of  manufacturing  variations. 

Although  binary  predictions  (shallow  or  deep)  are  convenient,  it  is  often  useful  to 
rank  injuries  on  a  continuous  scale,  so  that  the  size  of  the  burn  and  patient  comfort  often 
may  be  weighted  against  the  estimated  healing  potential.  The  proposed  algorithm  can 
be  used  to  produces  a  continuous'  scale  (O)  as  well  as  binary  predictions.  For  sites 
which  were  known  to  have  healed  in  the  third  and  fourth  weeks  post  burn,  there  is 
roughly  a  linear  correlation  between  instrument  response  and  healing  time  (Figure  40). 
The  estimated  number  of  days  required  for  healing  {H)  is 

H  =  30  0  +  15  (40) 

This  relationship  is  not  applicable  to  sites  which  heal  in  less  than  two  or  more  than  four 
weeks,  but  such  information  would  rarely  be  important  in  management  decisions. 

For  sites  which  heal  between  14  and  28  days,  the  error  in  this  model  (±  3  days) 
is  comparable  to  the  uncertainty  of  the  reference  method  (±  2  days).  Whether  or  not  a 
wound  was  healed  (entirely  reepethillalized)  is  difficult  to  determine  precisely.  Even 
among  the  Burn  Center  staff,  opinions  differed  by  a  day  or  two.  For  wounds  whose 
outcomes  were  described  by  patients  over  the  phone  or  extrapolated  from  the  nurses’ 
last  observations,  the  uncertainty  was  about  ±3  days. 

CONCLUSIONS 

With  the  LT  spectrometer  we  are  able  to  observe  absorbance  bands  from  most  of 
the  major  constituents  of  skin.  The  instrument  is  sensitive  enough  to  monitor  changes 
as  small  as  those  induced  by  only  a  few  degrees  elevation  in  surface  temperature.  The 
techniques  used  are  completely  non-invasive,  non-contacting,  and  non-destructive  and 
require  only  a  few  minutes,  making  them  suitable  for  use  on  burn  patients. 
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In  spite  of  the  non-linear  effects  of  light  scattering  and  multiple  layers,  hemoglobin 
speciation  can  be  estimated  from  visible  reflectance  measurements  of  burn  wounds. 
Further,  these  estimates  correlate  to  the  depth  of  the  injury.  Severely  elevated  levels  of 
methemoglobin  have  been  observed  in  deep  partial  thickness  and  full  thickness  burns. 
This  suggests  that  the  Imaging  Burn  Depth  Indicator  estimates  burn  depth  by  analyzing 
the  small  color  changes  caused  by  the  presence  of  methemoglobin.  The  IBDI,  initially 
designed  to  measure  deoxyhemoglobin,  can  be  significantly  improved  by  optimizing  its 
filters  to  detect  methemoglobin.  The  resulting  instrument  would  be  fast,  accurate,  non- 
invasive,  and  non-contacting  and  could  produce  a  map  of  a  wound,  color-coded  to 
indicate  healing  potential.  This  improved  instrument  could  be  used  as  a  practical  guide 
for  burn  treatment  decisions  and  may  assist  in  the  study  of  burn  physiology. 
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FUNDAMENTAL  STUDIES  OF  WATER 
INTRODUCTION 

In  order  to  understand  the  real  factors  that  lead  to  the  changes  in  reflection  spectra 
from  the  burn  sites,  we  have  to  investigate  the  temperature  dependence  of  the  water 
combination  bands  in  the  wavelength  range  used  in  our  in  vivo  studies  of  human  tissue. 
The  main  goal  of  these  studies  is  to  understand  the  structure  of  the  observed  peaks  and 
the  kinetics  of  the  real  molecules  and  structures  that  causes  the  observed  changes  in 
spectra.  Another  problem  arises  from  the  nonlinear  dependence  of  the  reflection  spectra 
on  the  extinction  and  scattering  of  the  tissue.  In  order  to  build  a  realistic  model  for 
temperature  dependent  reflection  spectra,  we  need  to  find  the  “pure  components"  of  the 
spectra  and  their  behavior  upon  change  of  temperature.  As  long  as  water  is  the 
dominant  component  of  tissue,  understanding  the  structure  of  water  and  specifically  how 
hydrogen  bonding  affects  NIR  spectra  is  one  of  the  prominent  problems  in  biophysics. 

The  structure  of  liquid  water  remains  a  most  puzzling  problem  in  the  chemistry  of 
solutions  and  electrolytes.  Despite  numerous  attempts  to  build  a  consistent  model  of 
water  through  empirical  data  collected  over  the  last  60  years,  it  must  be  admitted  that  no 
single  description  satisfactorily  explains  the  majority  of  the  collected  data.'*^  In  principle, 
vibrational  spectroscopy  provides  a  precise  tool  for  elucidation  of  the  various  species 
existing  in  water  and  water  solutions.  However,  even  in  such  a  straightforward  endeavor 
as  deriving  thermodynamic  properties  from  IR  absorption  spectrophotometry  and  Raman- 
spectroscopy,  the  results  are  in  poor  agreement  with  each  other  and  with  other  methods 
such  as  calorimetry.  The  main  source  of  this  ambiguity  in  spectra  interpretation  is 
concerned  with  the  problems  of  the  baseline  determination^®  and  with  the  difficulty  of 
interpolation  and  deconvolution  of  broad  bands  which  undergo  subtle  changes.  A 


j 


36 


number  of  investigators  have  attempted  to  deconvolve  the  absorption  bands  into  an 
arbitrary  number  of  gaussian  bands.  There  is  no  theoretical  reason  for  this  procedure 
and  the  results  are  not  consistent  with  any  known  theory.  Recent  advances  in  short 
wavelength  NIR  spectroscopy,  driven  by  new  technology  for  multichannel  silicon 
detectors,  makes  this  region  of  the  water  spectra,  where  only  overtone  and  combination 
bands  exist,  very  attractive  for  investigation.  Indeed,  many  models  of  water  explaining  dif¬ 
ferent  features  of  the  SW-NIR  exist  in  the  literature.'*®'®^ 

The  first  studies  of  NIR  combination  bands  were  done  more  than  three  decades 
ago  by  Suhrmann  and  Breyer.*®’^^  More  recent  investigation  of  NIR  bands  and  a 
suggestion  of  a  model  containing  three  water  species  have  been  described  by 
G.R.Choppin  and  K.Buijs  in  the  early  eO’s.^®"®®  They  proposed  that  three  bands  can  be 
distinguished  in  the  absorption  band  of  water  between  1100  and  1300  nm,  giving  the 
spectral  assignments  and  the  extinction  coefficients  for  each  band.  The  concentration  of 
each  of  the  three  absorbing  species  was  calculated  as  a  function  of  water  temperature. 
Using  this  model  and  the  concentration  of  each  species,  a  number  of  the  properties  of 
water  were  calculated.  The  species  were  assumed  to  be  non-bonded  (monomer)  and 
hydrogen  bonded  (dimer  and  trimeV)  intramolecular  complexes.  Their  enthalpies  were 
calculated  from  Van’t  Hoff  equilibrium  equations.  Subsequently,  different  authors 
proposed  models  of  water  with  two®®  and  even  five  species®®  and  many  theoretical 
papers  examining  different  models  appeared  in  the  last  three  decades.  The  most 
popular  is  the  model  proposed  by  Choppin  and  Buijs,  and  later  elaborated  by  Senior  and 
Vand*^,  that  explains  the  spectroscopic  data  in  terms  of  bonded  and  nonbonded  groups 
of  molecules  in  water,  is  referred  to  as  a  mixture  model.  Opposed  to  this  model  are 
continuous  models  of  the  water  structure.®®"®^  The  proponents  of  continuum  theories 
emphasize  the  continuous  evolution  of  the  broad  bands,  rather  than  trying  to  deconvolve 
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them  into  substructures.  From  our  point  of  view,  the  discrepancy  between  these  two 
classes  of  the  theoretical  models  appears  to  be  more  semantic,  and  they  are  supplemen¬ 
tary  rather  than  contradictory.  In  recent  study  of  infrared  spectra  of  aqueous  dispersions 
Hubner  et  al.“  have  pointed  out  the  danger  of  interpreting  complex  spectral  bands  solely 
on  the  basis  of  spectral  shifts.  They  cite  examples  where  deconvolution  of  the  complex 
band  explains  the  spectral  shift  in  terms  of  discrete  components,  while  the  frequency  shift 
of  overall  line  contour  leads  to  the  opposite  conclusions  about  the  phase  transition 
studied.  The  main  argument  against  the  structural  (mixture)  models  of  water  is  based 
on  a  statement  that  those  models  are  not  able  to  explain  the  entire  collection  of  empirical 
data  from  spectroscopy,  thermodynamic,  and  neutron  diffraction.  Nevertheless,  in  some 
recent  papers®®"®^,  spectroscopic  studies  of  the  2v,  +  Vg  combination  band  (960  nm)  in 
a  wide  range  of  pressures  and  temperatures  has  led  to  a  reconsideration  of  the  mixture 
model  involving  three  species.  In  these  papers,  certain  features  of  the  continuum  model 
are  retained  by  use  of  the  concept  of  energy  bands®^  rather  than  discrete  sharp  energy 
levels.  This  model  proposes  a  distribution  of  the  hydrogen  bond  energies  grouped 
around  three  component  bands  So.S^  and  Sg,  i.e.  a  model  featuring  a  relatively  small 
number  of  distinguishable  molecular  species  with  continuous  distribution  of  hydrogen 
bond  length  and  angles  associated  with  each  specie.  The  combination  band  peak  has 
been  resolved  into  three  gaussians  with  constant  positions  and  widths.®®'®’^  Similar  results 
have  been  presented  by  V.Fornes  gind  J.Chaussidon®®  for  the  Vj  +  V3  first  combination 
tone  of  water  molecules.  Both  groups  elaborate  quantitative  results  for  the  energy  of  the 
hydrogen  bond  formation  (rupture).  All  of  the  above  mentioned  investigators  derived  their 
data  by  fitting  the  spectral  bands  to  the  sum  of  tvyo  or  more  gaussians,  in  most  cases 
with  the  aid  of  an  analog  "curve  resolver". 

Recent  advances  in  multivariate  statistics  and  low-noise  NIR  spectrometers 
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provides  an  impetus  for  restudy  of  the  spectra  of  pure  water  as  a  function  of  temperature. 
Coupled  with  advances  in  instrumentation  for  high  precision  data  acquisition  are  advances 
in  data  analysis.  At  present,  there  seems  to  be  no  existing  methodology  for  quantitatively 
deciding  the  number  of  species  contributing  to  water  spectra.  Clearly,  arguments  based 
upon  the  existence  or  absence  of  isosbestic  point  do  not  appear  definitive  nor  do 
observations  of  spectral  shifts  of  the  mean  position  of  complex  peaks. 

It  is  the  purpose  of  this  study  to  describe  a  technique  for  obtaining  a  lower  bound 
to  the  number  of  species  contributing  to  a  series  of  spectra  taken  as  a  function  of  some 
external  variable(e.g.  temperature,  pressure)  within  the  linear  additivity  constraint.  This 
technique  which  we  call  "Chemical  Regression"  was  first  proposed  by  Box  and  later 
described  fully  by  Lawton  and  Sylvestre.®®’“  The  method  to  be  described  more  fully  in 
the  theory  section  relies  on  a  description  of  the  spectra  in  form  of  linear  combination  of 
the  eigenvectors  of  the  spectral  data  matrix.  Criteria  are  available  to  determine  a  lower 
limit  to  the  number  of  eigenspectra  required  to  describe  the  signal  variance  while 
eliminating  the  noise  variance.  If  the  number  of  eigenspectra  is  small  it  can  be  argued 
that  the  mixture  model  has  merit.  On  the  other  hand,  a  finding  that  a  large  number  of 
basis  vectors  must  be  retained  might  be  used  to  bolster  the  continuum  model.  Not  only 
can  a  lower  bound  to  the  number  of  components  be  found,  but  also  estimates  of  the 
spectra  associated  with  each  species  can  be  obtained.  This  is  accomplished  by  rotation 
of  the  abstract  eigenvectors  into  a  set  of  vectors  which  obey  physical  constrains 
appropriate  to  the  problem.  Thus  the  necessity  to  assume  that  the  spectral  shapes  of  the 
compo'^ents  is  gaussian  is  ei/minated  and  the  spectral  profiles  are  derived  directly  from 
the  data.  Constraints  may  be  imposed  not  only  on  the  spectra,  e.g.  positivity,  but  may 
also  be  imposed  upon  the  way  in  which  the  intensity  varies  with  the  external  variable  (i.e. 
we  may  postulate  a  model  for  the  effect  of  physical  variance  and  derive  the  physical 
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parameters  by  least  square  fitting).  Thus  it  would  appear  that  chemical  regression 
potentially  provides  further  insight  into  the  mechanism  for  spectral  changes  of  water  with 
temperature,  yielding  estimates  for  the  number  of  species  involved,  the  thermodynamic 
parameters  governing  their  interrelationship  and  estimates  of  the  spectra  of  each. 

In  this  study  we  have  used  this  new  method  to  reexamine  the  960  nm  combination 
band  in  pure  water.  This  band  demonstrates  a  strong  dependence  on  temperature 
changes  and  has  a  well  defined  isosbestic  point,  that  could  be  a  sign  of  the  presence  of 
two  or  more  species  in  pure  water.  Another  important  reason  for  choosing  the  SW-NIR 
range  for  examining  the  structural  model  of  water  is  connected  to  the  evident  fact  that  in 
this  range  the  shift  between  the  spectral  components-species,  if  they  do  exist,  would  be 
more  distinguishable^^,  and  their  separation  could  be  less  ambiguous. 

EXPERIMENTAL 

The  measurements  were  made  with  an  in-house  constructed  SW-NIR  diode-array 
spectrometer®’  and  a  4-cm  optical  path  length  quartz  cell,  and  repeated  with  the  Hewlett 
Packard  diode  array  spectrometer  HP-6582A  and  a  2-cm  cell  inside  the  HP-spectrometer. 
In  the  case  of  the  home-made  spectrometer,  the  radiation  was  delivered  to  the  cell  with 
a  fiber  optic,  which  made  temperature  control  more  convenient  and  reduced  fog  on  the 
windows.  Temperature  control  was  achieved  by  equilibrating  the  sample  cell  in  a  water 
bath  at  a  fixed  temperature.  The  temperature  of  the  cell  was  measured  with  a 
thermocouple  and  the  signal  was  digitized  and  registered  on  the  computer’s  hard  disk 
simultaneously  with  the  digitized  spectroscopic  data  from  NIR  spectrometer.  The 
measured  temperatures  are  believed  to  be  correct  to  within  0.5®C.  The  spectra  of  the 
pure  distilled  and  deionized  water  in  the  spectral  range  850-1 100  nm  were  taken  over  the 
temperature  range  10-80  °C.  The  cell  filled  with  CCI4  was  used  as  reference.  The 
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spectra  are  depicted  in  Figure  41 .  The  spectra  contain  visible  baseline  offset  and  a  slope. 
The  nature  of  this  bias  and  it’s  subtraction  are  discussed  below.  Nevertheless,  we  did 
not  take  special  measures  to  improve  the  absolute  measurements  of  the  temperature,  but 
the  rate  of  the  spectral  scanning  (60  scans  per  second)  allowed  us  to  reduce  the  noise 
by  repetitive  averaging.  We  believe  that  our  experimental  data  for  the  studied 
combination  band  are  the  most  precise  at  this  time. 

RESULTS 

The  temperature  dependent  offset  could  be  referred  to  the  changes  of  refractive 
index  of  water  with  the  temperature.  Fresnel  reflection  R  changes  according  to  the 
expression: 

R=(n-1)7(n+1)*  (41) 

where  n=nq/n^ ,  n^  and  n^  are  refractive  indexes  of  quartz  and  water  respectively.  The 
changes  in  n^  in  the  entire  temperature  range  are  less  than  1.5%  and  corresponding 
change  in  R  are  of  order  of  0.01 .  Taking  into  account  that  the  cell  has  two  interfaces  with 
water  and  the  effective  absorbance  will  increase  with  the  temperature  (the  refractive  index 
and  the  reflectance  are  respectively  decreasing),  the  baseline  shift  is  estimated  to  be 
within  0.02  absorbance  units,  in  good  agreement  with  observed  data.  Now, 
understanding  the  nature  of  the  baseline  offset  we  can  subtract  it.  The  bias  that  is  due 
to  the  strong  absorption  bands  existing  in  the  vicinity  of  the  investigated  peak  could  be 
approximated  by  the  straight  line  and  subtracted  in  the  same  way.  After  subtracting  the 
offset  and  the  slope  the  spectra  appear  more  regular  (Figure  42).  We  have  taken  into 
account  the  changes  of  the  water  density  with  the  temperature.  These  small  corrections 
are  within  1.5%  of  the  absorbance  value. 

In  order  to  build  a  realistic  model  of  water  we  need  to  determine  the  number  of 
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latent  variables  that  can  describe  the  variance  in  the  obtained  spectra.  There  are  different 
approaches  to  this  problem.  Most  of  them  are  based  on  the  mathematical  procedure  of 
separating  the  noise  or  the  bias  from  the  real  variance.  This  mathematical  procedure  is 
called  Singular  Value  Decomposition  (SVD).  It  consists  of  representation  of  the  spectral 
matrix  of  absorption  A  in  the  form: 

A=USV’  (42) 

where  S  is  diagonal  matrix  with  the  nonnegative  diagonal  elements  in  decreasing  order, 
and  U  and  V  are  unitary  matrices,  V’  denotes  the  transposed  V  matrix.  The  columns  of 
A  represents  the  spectral  vectors  at  a  fixed  temperature.  The  rows  of  A  contain 
information  about  the  variation  in  absorption  at  a  fixed  wavelength.  We  can  also 
represent  the  experimental  matrix  A  in  the  form 

A=DC  +  E  (43) 

where  matrix  D  contains  only  spectral  components  of  the  mixture  and,  the  matrix  C 

represents  the  concentration  of  the  different  components  depending  on  the  variable 
parameter  (in  our  case  the  temperature).  Matrix  E  contains  experimental  noise.  From 

Equations  42  and  43  we  can  derive  an  equation: 

DC  =  USV’  (44) 

where  U,  S,  and  V  are  truncated  matrices  containing  only  the  vectors  and  eigenvalues 
which  are  meaningful  for  the  problem,  namely  they  take  into  account  the  variance  in  the 
spectral  matrix  A  and  reject  the  noise.  The  representation  of  the  experimental  matrix  A 
in  the  form  of  Equations  43  and  44  is  well  known  as  Principal  Component  Analysis 
(PCA).“  PCA  is  a  powerful  tool  in  separating  the  noise  from  the  real  spectral 
components  and  estimating  the  most  likely  number  of  latent  variables  (meaningful  singular 
eigenvalues)  in  the  matrix  A,  that  can  describe  the  variance  in  A  with  reasonable 
precision.  Since  the  determination  of  the  number  of  "real"  Principal  Components  and 
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elimination  of  those  which  describe  the  variance  of  the  noise  is  an  ambiguous  task,  we 
use  auxiliary  techniques,  such  as  Evolving  Factor  Analysis®^.  The  main  concept  of  EFA 
is  a  graphical  representation  in  which  normalized  (divided  by  the  trace)  singular  values 
of  the  submatrices  of  matrix  S  in  descending  order  versus  the  constituent,  are  displayed 
and  examined  for  "emerging"  eigenvalues.  EFA  has  been  used  to  analyze  spectra, 
consisting  of  strongly  overlapped  peaks®^.  The  normalized  eigenvalues  usually  are 
depicted  on  a  logarithmic  scale.  For  the  matrix  A,  corresponding  to  the  spectra  shown 
in  Figure  42,  EFA  evidently  shows  three  emerging  eigenvalues.  The  results  of  EFA  are 
illustrated  in  Figure  43.  The  investigation  of  the  next  spectral  eigenvectors  from  matrix  U 
shows  that  they  contain  much  more  noise  than  the  first  three,  and  hence  they  are  useless 
for  constructing  the  set  of  real  spectral  components.  Therefore,  PCA  analysis  suggests 
the  same  rank  of  three.  Consistent  with  the  estimation  of  rank  three  is  the  structure  of 
the  spectral  eigenvectors  as  shown  in  Figure  44.  These  eigenvectors  are  largely  free  of 
noise.  Examination  of  the  higher  eigenvectors  showed  them  to  be  highly  contaminated 
with  noise  (Figure  45)  and  too  irreproducible  to  be  of  value  in  representation  of  our  data 
set.  Hence  within  our  current  signal  level  and  reproducibility,  a  basis  set  of  three 
components  seems  defensible.  Now,  starting  with  the  assumption  of  three  components, 
we  can  try  to  rotate  our  spectral  eigenvectors  (Figure  44)  to  the  set  of  real  spectral  peaks. 
The  procedure  of  nonlinear  regression  based  on  the  postulated  reaction  has  been 
described  by  Shrager®^  ®®.  This  procedure  of  recovering  the  real  spectral  components 
and  concentrations  (titration),  namely  calculating  the  matrices  D  and  C  from  the  matrices 
U,S  and  V  implies  that  one  makes  some  reasonable  assumptions  about  the  kinetics  of 
the  system  under  consideration,  and  then  finds  the  transformation  which  fits  the  rotated 
vectors  of  the  V-matrix  (so  called  "scores")  to  the  imposed  model.  After  making  the 
decision  about  the  rank  of  the  data,  the  next  step  is  to  postulate  a  model  for  the  kinetics 


of  the  reaction.  This  will,  in  turn,  provide  the  means  for  recovery  of  the  pure  spectral 
components  and  the  equilibrium  constants.  From  the  above  described  analysis  of  the 
rank  of  our  spectroscopic  data,  it  follows  that  there  are  three  distinguishable  molecular 
species  in  pure  water.  The  postulation  of  thermodynamic  equilibrium  among  all  three 
species  leads  to  a  set  of  equations  for  the  concentrations  Cj  of  different  species  : 

C2  =  C,exp(-AHi2/RT  +  AS^ 

C3  =  C^exp(-AHi3/RT  +AS3)  (45) 

+  C2  +  C3  =  1 

where  the  AHy  and  ASj  are  the  relevant  enthalpies  and  entropies  of  the  imposed  reaction. 
The  equilibrium  constants  are:  Ki=C2/C^,  K2=Ci/C3  ,  K3= 03/03= K,K2  ,  and  log(Kj) 
versus  reciprocal  temperatures  1/T  are  linear  functions.  Under  these  assumptions  we 
can  rotate  the  titration  eigenvectors' from  V  with  the  aid  of  3  by  3  matrix  T,  in  order  to 
yield  the  equation: 

V’=TO  (46) 

Comparing  Equations  42  and  44  we  can  derive  the  matrix  D,  containing  the  pure  spectral 
components, 

D=UST  (47) 

The  problem  of  calculating  the  transformation  matrix  T  is  central  to  this  nonlinear 

regression  procedure.  We  used  the  Nelder-Mead  simplex  optimization  algorithm,  that 
solved  Equation  46  with  four  unknown  nonlinear  parameters  A  H,2,  A  H13,  ASg  and  AS3,  and 
nine  unknown  linear  parameters,  that  actually  determine  the  matrix  T,  in  the  least  square 
sense.  The  routine  returns  the  unknown  parameters  and  the  norm  of  the  difference 
vector  between  the  functions  Cj  from  Equation  45  and  rotated  vectors  V,  that  serves  as 
a  measure  for  this  nonlinear  fitting.  The  original  "scores"  from  the  V  matrix  and  rotated 
to  the  real  concentration  Cj  are  depicted  in  Figures  46  and  47  respectively.  Now  after 
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calculating  the  T  matrix  we  can  recover  the  pure  spectral  components  with  the  aid  of 
Equation  47.  Three  spectral  components  shown  in  the  Figure  48  were  recovered  with  the 
above  described  nonlinear  regression  technique.  Figure  49  shows  the  spectral 
components  recovered  from  different  data,  obtained  from  different  spectrometers,  and 
demonstrates  the  level  of  confidence  in  recovering  the  spectral  components  with  the 
chemical  regression.  The  returned  enthalpies  for  the  reactions  AHjj  are:  AHi2=2.6 
kcal/mol,  =3.6  kcal/mol.  The  first  one  corresponds  to  the  well  known  value  nf  the 
heat  of  hydrogen  bond  formation  in  pure  water.  The  latter  is  probably  the  heat  of 
formation  of  the  intermediate  form  of  hydrogen  bonded  specie  in  water.  This  dependence 
is  very  sensitive  to  the  recovered  nonlinear  parameters  and  should  viewed  qualitatively 
rather  than  reflecting  real  entropies  and  enthalpies  of  the  hydrogen  bond  formation. 
Stability  analysis  of  this  nonlinear  regression  procedure  have  been  made,  by  adding 
artificial  random  noise  to  the  experimental  spectral  data.  The  routine  of  chemical 
regression  demonstrates  high  stability  in  recovering  spectral  components  but  fairly  poor 
reproducibility  for  concentrations.  This  feature  of  chemical  regression  is  not  unknown  but 
has  to  be  understand  in  terms  of  illposed  problems®®.  The  values  cited  above  correspond 
to  the  spectral  components  depicted  in  Figure  48  and  are  probably  not  very  reliable,  but 
the  shape  of  the  recovered  spectral  components  are  reproducible  for  different  sets  of 
data  and  demonstrates  high  stability  upon  small  random  perturbation  of  the  initial  spectral 
data.  Using  pseudoinverse  operator  and  recovered  spectral  components  we  were  able 
to  calculate  the  experimental  concentration  matrix  (closed  circles  in  Figure  47)  from 
the  matrix  D. 

=pinv(D)A  (48) 

pinv  in  Equation  48  denotes  operation  of  pseudoinversion,  according  to  definition 
pinv(A)A  =  i,  were  I  is  a  unity  matrix. 


The  spectrum  which  belongs  to  the  specie  which  dominates  at  high  temperature 
is  the  most  blue  shifted  and  narrowest  of  the  spectral  components,  indicating  the  least 
amount  of  hydrogen  bonding.  The  spectrum  from  the  specie  which  exist  through  the 
range  of  temperatures  is  intermediate  in  spectral  shift  and  width.  The  low  temperature 
specie  is  the  most  red  shifted  and  broadest  consistent  with  a  high  degree  of  hydrogen 
bonding,  and  reflects  the  statistical  nature  of  the  different  length  of  hydrogen  bonds  in 
water  clusters.  Thus  the  evolution  of  the  spectra  is  in  agreement  with  the  many  studies 
of  effect  of  hydrogen  bonding  on  spectra. 

Our  analysis  of  the  temperature  dependent  spectra  of  the  960  nm  combination 
tone  peak  demonstrates  that  the  three  spectrally  and  structurally  distinguishable  species 
are  needed  to  explain  the  evolution  of  the  spectra  within  the  framework  of  a  simple 
thermodynamic  model.  The  attempt  to  deconvolve  the  spectra  imposing  a  two  species 
model  (bonded  and  nonbonded  structures,  being  in  thermodynamic  equilibrium)  derives 
unrealistic  spectra  of  the  components  and  thermodynamic  parameters.  It  seems  very 
probable  that  more  precise  spectral  data  and  absolute  temperature  measurements  would 
lead  us  to  reconsideration  of  the  minimal  number  of  species  needed  and  consequently 
the  imposed  thermodynamic  model.  Such  studies  would  demand  a  substantial 
improvement  of  the  spectroscopic  apparatus  and  probably  a  more  diverse  approach  to 
this  problem,  exploiting  different  methods  of  CARS  spectroscopy  of  Raman-active 
vibrational  transitions  or  more  precise  investigation  of  lower  overtones  in  the  NIR  range 
of  water  spectra.  Such  studies  are  complicated  by  the  high  level  of  noise  due  to  the  high 
extinction  coefficient  and  problems  of  temperature  and  flow  control  inside  very  thin  cells. 

CONCLUSIONS 

In  this  study  we,  for  the  first  time,  applied  the  powerful  technique  of  multivariate 


statistics  and  nonlinear  regression  to  the  problem  of  water  structure.  Our  findings  support 
the  mixture  model  of  water  with  at  least  three  species.  For  the  first  time,  we  recovered 
the  pure  spectra  of  these  species  from  an  imposed  explicit  physical  model. 
Nevertheless,  at  this  time  we  don’t  have  a  clear  understanding  of  the  physical  nature  of 
these  species  and  the  relevant  energies  of  reactions.  However,  this  model  can  be  used 
as  an  analytical  approach  to  the  interpretation  of  water  spectra  and  even  for  measuring 
the  temperature  of  the  water-containing  samples.  The  obtained  spectral  components  can 
be  exploited  for  developing  a  statistical  theory  of  hydrogen  bonding  in  water  and  can  be 
helpful  in  verification  of  the  numerous  theories®®  ®^  which  try  to  calculate  the  shape  of  the 
absorption  and  Raman  band  in  pure  water.  The  pure  components  found  from  the 
temperature  dependent  spectra  of  pure  water  could  be  used  for  interpretation  of  the 
temperature  dependence  of  the  NIR  spectra  of  human  tissue  and  for  improving  of  the 
quantitative  prediction  of  the  depth  of  burn  wounds. 
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MONTE  CARLO  SIMULATION  OF  HETEROGENEOUS  TISSUE 

The  study  of  light  propagation  and  reflection  from  tissue  is  central  to  many  medical 
and  biomedical  applications  of  light.  For  diagnostic  purposes,  such  as  noninvasive  optical 
spectroscopy,  the  light  that  is  diffusely  reflected  from  or  transmitted  through  the  tissue 
may  be  measured  to  probe  the  metabolic,  physiologic,  or  possibly  the  structural  status 
of  the  tissue.  In  previous  sections  of  this  report,  we  stated  that  Beer’s  laws  and 
multivariate  were  used  to  analyze  reflection  of  light  from  tissue,  a  medium  in  which 
scattering  and  absorption  occur.  Beer's  Law  is  particularly  troublesome  when  both 
scattering  and  absorption  are  strong.  Any  further  improvement  of  our  predictions  based 
on  a  multivariate  statistics  will  depend  on  a  proper  algorithm  of  reflectance  spectra 
evaluation,  involving  subtraction  of  temperature  dependent  water  spectra  and  comparison 
of  the  transmittance  in  vivo  and  reflectance  spectra  of  oxy  -  and  deoxy-  hemoglobin. 

The  problem  of  the  deconvolution  of  the  remittance  spectra  of  tissue  have  been 
discussed  in  many  papers.®®  It  is  well  known  that  changes  in  the  distribution  of  scattered 
light  in  the  skin  that  result  from  variations  in  the  hematocrit  and  volume  fraction  of  blood 
contained  in  the  dermis  may  affect  the  calibration  of  the  instruments  used  in  pulse  blood 
oximetry  or  burn  wound  evaluation.  It  could  be  extremely  difficult  to  perform  controlled 
experimental  studies  of  these  effects,  especially  in  the  frequency  domain.  Therefore 
mathematical  models  and  computer  simulations  of  light  propagation  and  diffusive 
reflection  from  tissue  are  essential  to  guide  current  applications  and  to  promote  further 
developments  in  biomedical  optics. 

The  three  most  common  mathematical  approaches  to  the  study  of  light 
propagation  in  tissue  are  based  on  the  Kubelka-Munk,  random  walk,  and 
radiative-transport  model.  Most  of  the  studies  have  been  done  for  the  semi-infinite 
homogeneous  model  of  tissue,  and  for  fixed  scattering  parameters  and  a  given 
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wavelength  of  light.  The  importance  of  modeling  in  the  wide  wavelength  domain  follows 
from  our  simple  simulations  made  for  continuous  homogeneous  medium  with  the  aid  of 
analytical  expressions  obtained  by  Wilson.®®  We  have  used  the  spectra  of  oxy-,  deoxy-, 
and  met-  hemoglobin  obtained  during  our  in  vitro  studies,  and  have  calculated  the 
extinction  spectra  that  would  be  measured  by  our  methods  in  tissue.  The  original  and 
simulated  spectra  are  shown  in 

Figure  50.  Solid  lines  represent  original  spectra,  and  dashed  lines  represent  simulated 
tissue  spectra. 

Even  this  simplistic  approach  shows  the  significant  difference  in  spectral  shape, 
and  emphasizes  the  importance  of  evaluating  reflectance  spectra  properly  in  order  to 
improve  our  quantitative  interpretation. 

Our  plan  for  research  included  following  steps; 

1 .  Development  of  a  working  program  for  simulations  of  photon  transport  in  a  random 
scattering  heterogeneous  medium  with  inclusions  (i.e.,  veins)  with  different  scattering  and 
absorption  parameters. 

2.  Numerical  simulations  of  light  distribution  inside  and  outside  homogeneous  scattering 
medium  for  a  fixed  wavelength  and  for  a  wide  range  of  wavelengths. 

3.  Comparison  of  our  numerical  simulations  with  the  existing  experimental  and  theoretical 
studies. 

4.  Numerical  simulations  of  the  reflected  light  distribution  from  the  multilayer 
heterogeneous  medium  with  inclusions  for  the  wide  wavelength  range  of  visible  and 
SW-NIR  light  (so-called  "therapeutic  window"). 

5.  Numerical  and  theoretical  studies  of  the  pulse  response  from  the  heterogeneous 
media. 

We  were  not  able  to  fulfil  this  program  because  of  the  termination  of  funding,  but 


even  the  preliminary  results  appear  promising  and  informative. 

The  basic  principles  of  four  different  types  of  remittance  spectrophotometry  could 
be  understood  from  Figure  51.“  We  developed  a  working  program  in  FORTRAN 
language,  that  was  able  to  run  on  different  computers;  VAX/VMS,  IBM-6000/RS,  and 
HP-6000.  The  program  was  able  to  simulate  the  random  walk  of  photons  in  a  scattering 
medium  with  inclusions  taking  into  account  such  effects  as  internal  reflection  and 
refraction  due  to  the  different  indices  of  refraction  of  tissue  and  air.  Monte  Carlo 
simulations  that  include  effects  of  geometrical  optics  have  not  yet  been  discussed  in  any 
of  the  numerous  published  studies.  We  used  as  a  guideline  for  our  Monte  Carlo 
simulations  a  model  by  Bonner^o.  Basic  steps  in  the  calculation  of  the  single  photon 
path  are  obvious  from  Rgure  52.  The  term  "weighting"  in  Figure  52  refers  to  the 
"intensity"  of  the  photons,  which  decreases  with  path  travelled  according  to  the  assumed 
absorption  coefficients  in  the  medium.  The  computer  code  in  FCRTRAN  is  listed  in  the 
Appendix  B.  The  program  was  structured  in  several  subroutines  that  made  it  flexible  for 
fast  changes  of  parameters,  shape  of  the  light  source,  and  the  number  and  shape  of 
inclusions.  The  shape  of  inclusions  could  be  given  by  algebraic  equation  or  numerically. 
We  used  basically  two  types  of  light  sources,  point  and  a  radial  source  with  a 
homogeneous  distribution  of  intensity.  The  typical  number  of  inserted  photons  that  gave 
a  reproducible  intensity  distribution  inside  and  outside  the  tissue  was  100000.  The 
simulations  were  time  consuming.  Cne  run  (for  the  fixed  set  of  scattering  parameters) 
took  typically  6-8  hrs.  on  the  VAX  computer  and  0.5  hours  of  CPU  time  on  an  IBM-6000 
workstation  with  performance  of  2-3  Mflops  per  second. 

The  model  of  the  heterogeneous  tissue  with  radial  inclusions  is  shown  in  Figure  53. 
The  results  of  the  simulations  were  plotted  as  intensity  contour  lines.  Figures  54-56 
demonstrate  the  intensity  distribution  for  a  point  light  source  just  inside  the  homogeneous 
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medium  for  the  case  without  inclusions,  one  radial  inclusion,  and  two  radial  inclusions, 
respectively.  The  scattering  parameters  are  shown  on  the  plots.  The  intensity  distributions 
for  a  radial  source  in  the  homogeneous  medium  without  inclusions  and  with  one  radial 
inclusion  are  depicted  in  Figures  57-58.  Figures  59-60  demonstrate  the  spatial  and 
temporal  behavior  of  the  reflected  light,  corresponding  to  Figure  57. 

Therefore,  in  this  chapter  we  report  an  innovative  computer  code  for  Monte  Carlo 
simulations  of  light  remittance  from  heterogeneous  media.  Our  results  obtained  over  a 
comparably  short  period  of  time  with  comparably  modest  computer  resources  could  be 
an  impetus  for  further  studies  in  this  new  and  promising  field  of  biomedical  research. 
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TABLE  X 


Table  of  Assignments. 


1764  Fat 
1731  Fat 
1728  Collagen 
1708  Fat 
1689  Collagen 

1460  Water  (0th  derivative) 
1400  Water  (2nd  derivative) 
1275  Collagen? 

1213  Fat 
1189  Fat 

1188  Protein  in  muscle  tissue 
1153  Water 
1038  Fat 
1014  Fat 
958  Water 
929  Fat 

760  Deoxyhemoglobin 
576  Oxyhemoglobin 
555  Deoxyhemoglobin 
541  Oxyhemoglobin 


CHj  stretch,  1st  overtone 
CHj  stretch,  1st  overtone 

CHj  stretch,  1st  overtone 

OH  stretch,  1st  overtone 


CHj  stretch,  2nd  overtone 
CHj  stretch,  2nd  overtone 

OH  combination 
CHj  combination 
CHj  combination 
OH  stretch,  2nd  overtone 
CHj  stretch,  3rd  overtone 
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Rgure  la:  Visible  Spectra  of  Hemoglobin  Derivatives 


.  . 

Cjyf  boxy  W  « loWt^ 
—  *  —  *  ^y£i.w.ome+i\#,»\r^lcb/r 


Hb02 

Hb 


Rgure  1  b:  Short-wave  Near-infrared  Spectra  of  Hemogiobin 
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Figure  23.  Loadings  from  principal  component  analysis  of  baseline-subtracted 
spectra.  Eigenvectors  a)  1  and  2,  b)  3  and  4.  Eigenvectors  (-)  1  and  3  and  (...) 
2  and  4. 
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Figure  28.  Scores  from  principal  component  analysis  of  normalized  spectra,  (x) 
shallow  and  (o)  deep  injuries. 
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Hemoglobin  Spectra  Model,  fitted  spectra 
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Figure  37.a  Burn  Depth  Indicator  Simulation.  (— )  a  typical  spectrum  of  a  deep 
burn,  and  (•  •  •)  the  LED  output  functions. 
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Figure  37.b  Imaging  Burn  Depth  Indicator  Simulation.  (— )  a  typical  spectrum  of 
a  deep  burn,  and  (•  •  •)  the  filter  transmittance  functions. 
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Rgure  38.  Simulated  instrument  responses,  (a)  the  burn  depth  indicator,  (b)  the 
imaging  burn  depth  indicator,  (x)  are  responses  for  shallow  burns,  (o)  for  deep 
burns.  (-)  is  the  best  discriminant  line  found  by  vanUew. 
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Fig.jT/  Types  of  remittance  measurements.  The  four  types  of  remittance 
measurements  are  1 )  steady-state  measurement  of  the  total  remittance  as 
either  reflectance  R  and/or  transmittance  T,  2)  time-resolved  measure- 
ment  of  total  remittance  /?(/ )  as  a  function  of  time,  3)  steady-state  mea¬ 
surement  of  the  local  spatial  distribution  of  remittance  R(  p )  along  the 
tissue  surface,  and  4)  the  time-resolved  measurement  of  the  local  spatial 
distribution  of  remittance  /t(p,  r)  along  the  tissue  surface  as  a  function 
of  time.  The  photon  pathlength  in  the  tissue  is  L  *  ct/n  where  c  is  the 
in  vacuo  speed  of  light  and  n  is  the  tissue  refractive  index. 


INTENSITY  (ARB.  UNITS) 


Appendix  A 

MATLAB  Code  for  Rotation  and  Single-Plane  Discrimination 


function  [Xmean, Scale, Rota]  =  spinc(X,C.Options) 

%  X  =  ail  data  (samples  by  variables,  more  samples  than  variables) 
%  C  =  class  vector,  containing  I’s  and  2's  only,  O’s  for  ignored  data 
%OUTUNE 
%  mean  center 
%  scale 
%  rotate 

%SUBROUTEINS  NEEDED 
%  std,  varwt,  fisherwt 
%OPnONS 

%  1  scale  (Range,  Auto,  Variance,  Fisher,  none) 

%  2  rotation  steps 

%  3  number  of  rented  variables  to  calculate 
%function  [Xmean,Scale,Rota]  =  spinc(X,C,Options) 


^  fO  /O  TO  TO  TO  TO  TO  TO  TO  TO  ^3  TO  Al  TO  TO  TO  TO  TO  TO  /O  TO  ^  ^  ^9  fO  ^3  rv  ^  n  /Q/O  rO  /O  ^  ^/OfO  ^yO  TO  A9 /O  TO  AO  /D  /O  AO  AO /D  AOTD  aO  ADTD 

SETUP 

f9  TO  TO  ^O  TO  TO  ^0  ^O  ^0  TO  f9  ^0  AT  TO  TO  ^O  'O  TO  ^0  TO  TO  TO  TO  fv  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  fv  TO  TO  AO  TO  AO  TO  TO  AO  A3  TO  A3  TO  AO  AO  TO  70  AO  AO 

^****************************'«***************«**************  qheck  sizes  ***** 

%CHECK  SIZES 
(m,n|  =  size(C); 

if  m-  *  1  &  n~  =  1,  error(‘Rotation;  Too  many  constituents.’),  end 
ifn>m,  C=C’:  Im,nJ=si2e(C):  end 
[np,nv|  =  size(X); 

lfnp-=m&nv-=m.  error('Rotatlon;  Incompatabie matrices’),  end 
ifnv==m,  X=X’:  (np,nv|=size(X);  end 
clear  m  n 

if  any(C~  =  1  &  C~  =2  &  C~  =0),  error(’Rotation:  Bad  constituent  values.’),  end 
n1  =  firKl(C=  =  1); 

if  length(nl)<2,  error(’Rotation:  Too  few  of  class  #1’),  end 
n2  =  find(C=  =2); 

if  Iength(n2)<2,  eiTor(’Rotation:  Too  few  of  class  #2’),  end 

if  ~any((1  2  3  4  5]  =  =  Options(1)), 
error(’Rotation;  Not  a  valid  scaling  option.’) 
end 

if  length(Options)<3,  0ptions(3)  =  nv-1;  end 
if  Options(3)  >  nv-1 ,  Options(3)  =  nv- 1 ;  end 


%SETUP  MATRICIES 


PhiStepSize  =  1  /  0ptions(2); 

Rota  =  eye(nv);  %Main  rotation  matrix 

Cyde  =  0; 


SETUP  MATRICES 


***** 


TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  ^  TO  TO  TO  TO  TO  TO  TO  TO  TO  TO  AO  TO  TO  TO  TO  TO  TO  TO  TO  TO  AOTOTOAOTOTOAOAO  TO  ^  TO  TO  TO  TO 

%  PREPROCESS  DATA 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%MEAN-CENTER 


Xmean  =  mean(X); 

Data  =  X  -  ones(np,l)*Xmean; 


SCALE 


%SCALE  (Sharaf.lllman, Kowalski:  "Chemometrics",  p.193) 


if  Options(l)  =  =  1.  % 

Range  scaling 

for  1=1  :nv,  Scale(i)  =  1  /  max(abs(Data(:.i)));  end 
elseif  Options(l)  =  =2.  % 

Autoscaiing 

for  i=l:nv,  Scale(i)  =  1  /  std(Data(;,i)):  end 
elseif  Options(l)  =  =3,  % 

Variance  weights 

for  i=1:nv,  DataTmp(:,i)  =  Data(:,i)/max(abs(Data(:,i)));  end 
[Scale]  =  varwt(DataTmp,C); 
fori=1;nv.  Scale(i)=Scale(i)/max(abs(Data{:,i))):  end 
elseif  Options(l)  =  =4,  % 

Fisher  weights 

for  i=l:nv,  OataTmp(:,l)=Oata(:,i)/max(abs(Data(:,i))),  end 
[Scale]  =  fisherwt(OataTmp,C); 
for  1=1  :nv,  Scale(i)=Scale(i)/max(abs(Oata(:.i)));  end 
elseif  Options(l)  =  =  5. 

%no  scaling 
end 


forl=1:nv,  Data(;,i)=Data(:,i)*Scale(i):  end 


/D  Vv  VO  /O  vO  TP  70  /O  TO  rO  70  TO  TO  TO  TO  TO  rO  fO  TO  /O/O^^Al/O^AO/OAO/OAO/O  AQTD  AO  TO  AO  TO  AOAO/O/OAOAOAOAO  AO  TO  AO  AO  AO  A>  AO  TO  AO 

%  ROTATE  (two  vectors  at  a  time) 

o/oz^Q4,<^c^c^c^OAc^Q^q/(^oAoz.Q;,q^Q^qL 

TO  AO  TO  AO  AO  /O  AO  ^  TO  /O  /O  AO  VO  ^  ^  /O  TO  AO  TO  AOAO^AOAOAOAOAO/O  AO  AOAOAOAOAO^AOAOAOAOAOAOAOAOAO/O  TO  /OAOAOAOAOAOAOAO/O 


for  van  =  1  ;Options(3). 
forvar2  =  (var1  +  1):nv, 


TWO-VECTOR  SUB-MATRICES 


%Choose  vectors 
Data2  =  Data(;,[var1,var2]); 

Rota2  =  Rota(;,{var1,var2i); 

Axis  =  sqrt(2)  *  (-1  1  -1  1]  *  max(max((Data2))); 


ROTATE  (loop) 


***** 


%Spin  about  [0  0] 
for  phlstep=  1  :Options(2), 

%+  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  -(-  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  -n-  +  +  +  +  -(-+  Rotation  Angle 
%Rotate  data 

phi  =  phistep  *  PhiStepSize  *  (2* pi); 
rot  =  [cos(phi),  -sin(phi);  sin(phi),  cos(phi)J; 
data  =  0ata2  *  rot; 
rota  =  Rota2  *  rot; 

%+  +  +  +  +  +  +  +  +  +  -I-  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  Distance  Between  Classes 
%Distance  between  classes 
[d.signj  =  vanNt(data(:,l),C); 

If  sign  =  =  d  =  -d;  end 
dist(phistep)  =  d; 
end  %for  phistep 


FIND  MAXIMUM  ANGLE 


***** 


)to 


%Maximize  distance 

%+  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  Angle 
maxphi  =  find(  dist  =  =  max(dist) ); 

if  length(nftaxphi)  >  1,  maxphi  =  maxphi(l);  end 
phi  =  maxphi  *  PhiStepSize  *  (2*pi); 
rot  =  [cos(phi),  -sin(phi);  sin(phi),  cos(phi)]; 

^****************************************«*********  SUB-MATRICES  ***** 

%New  sub-matrices 

%+  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  +  New  Matrices  +  +  +  +  + 
data  =  0ata2  *  rot; 
rota  =  Rota2  *  rot; 

Data(:,[var1,var2])  =  data; 

Rota(;,(var1,var2j)  =  rota; 

disp([’Rotation:  Rnnished  variables  ',int2str(varl),’  and  ’,int2str(var2)]) 
end  %for  var2  =  (van  +  1):nv 
end%forvar1  =  1:(nv-1) 


function  [DataNew]  =  spinpp(mean, Scale, Rota.Xp) 
%function  [DataNew]  =  spinp(Xmean.Scale.Rota,Xp) 

(np  ,nv]  =  sizepCp); 

Data  =  Xp  -  ones(np,l)*Xmean: 
fori=1:nv.  Data(:,i)=Data(:.i)*Scale(i):  end 
Data  =  Data  *  Rota; 

fori=1:nv,  Data(:,i)=Data(:.i)/Scale(i):  end 
DataNew  =  Data  +  ones(np,l)*Xmean: 


function  [discPoint]  =  disc(D,C); 

%  DISCc.m 

%  Finds  the  most  discriminating  point  to  separate  two  classes 
%  M.R.Balkenhol,  1/12/92 
% 

%  D  =  data  in  two  columns,  column  1  is  for  discriminating 
%  C  =  class  vector,  1  's  and  2's  only 
% 

%  minE  =  best  point  (index  and  data  value) 

%  E  =  number  of  misclassifications  [index,  ciassi,  class2,  both  classes] 
% 

%function  [discPoint]  =  disc(D,C); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
INPUT  and  SETUP 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
mInD  =  mln(D(:,1)): 
maxO  =  max(D(;,1)): 
maxi  =  300; 

Error  =  999*ones(maxl+1,l): 

Error2  =  999*ones(maxi+1,1); 

Index  =  2eros(maxi+1,1); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  LOCATION  OF  MOST  DISCRIMINATING  POINT 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

^****«**********«******«*********«************«***  Mispredictions  ***** 

for  i=1;maxi, 

I  =  i*(maxD-minD)/rnaxi  +  minD; 
el  =  sum(  D(:,1)>=j  &  C=  =  1  ): 
e2  =  sum(  D(:,1)<  j  &  C==2  ); 

Error(i)  =  el  +  e2; 
end 

[minEs]  =  min(Error); 
minis  =  flnd(Error=  =  minEs); 

()^*********************************«*******  Smallest  Sum-of-Squares  Error  ***** 

l=mln(minls) :  max(minls); 

fori=l, 

j  =  i*(maxO-minO)/maxi  +  minD; 
n1  =  flnd(D(:,1)>=j  &C=  =  1  ); 
n2  =  flnd(  D(:,1)<  j  &C==2); 

Error2(i)  =  sum(  ( D([n1;n2],1)  -j).'^  ); 
end 

[minE,minl]  =  min(Error2); 
n  =  length(mlnE); 

if  n  >  1,  minE  =  minE(  round(n/2),  : );  end 


discPoint  =  minl*(maxD-minD)/maxi  +  minD; 


1^-3 


function  [Cest,ep]  =  discp(disc,Dp,Cp); 

%  DISCp.m 

%  Sorts  a  predictions  set  (Ciassi  <  disc,  Ciass2  >  disc) 

%  M.R.Baikenhol,  1/12/92 
% 

%  Dp  =  data  in  two  coiumns,  coiumn  1  is  for  discriminating 
%  Cp  =  ciass  vector,  I’s  and  2's  oniy  (optionai) 

% 

%  disc  =  discriminating  point,  minE(1,2)  from  disc.m 
% 

%function  [Cest]  -  discp(disc,Dp,Cp); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
iNPUT  and  SETUP 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
[np.nv]  =  size(Op); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  SORTiNG  DATA  POINTS  with  respect  to  DISC 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
for  i=1:np, 
if  Dp(i,1)  <  disc, 

Cest(i)  =  1; 
eiseif  Dp(i,1)  >  disc, 

Cest(i)  -  2; 
else 

Cest(l)  «=  0; 
end 
end 

Cest  =  Cest’; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  DISPLAY  (if  Cp  exists) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
if  nargin  >  2, 

errors  =  find(  Cest  ~  =  Cp  ); 
plot(find(Cest==1),  Dp(Cest=  =  l,1),  ’c6+’,.. 
find(Cest==2),  Dp(Cest==2,l), 

[1  npl,  [disc  disc],  ’b:’,.. 
errors,  Dp(errors,l),  ’c5o’) 
ep  =  length(errors); 
end 


PROGRAM  SVEIN 


yC  13 


REALNO^T 

DIMENSION  A(2),  FLUX(100,100)^LR(11,600)ALT{11,600) 

INTEOER*41X,IY 

OPEN(UNlT=l.nLE=’FVElN.C}GP^TATUS-’NEW) 

OPEN(UNIT=2JTLE=’LVEIN.GGP,STATUS-'NEW’) 

C  DEFDOTION  OF  PARAM  ATERSS 

C  FLUX  ARRAY  THAT  IS  CELL  HISTORIES  ARE  STORED  TO  CALCULATE 
C  A  GRAPHIC  DISPLAY  OF  FLUX. 

CAULALT  ARRAYS  THAT  RECORD  PHOTON  PATHLENGTH  ABOVE  AND  BELOW 
C  THE  TISSUE  INTERFACE 

C  THETA  ANGLE  THAT  PHOTON  WILL  MOVE  NEXT 
C  PLNGTH  PATH  LENGTH  THAT  PHOTON  TRAVELS  IN  ONE  STEP 
C  RHO  DENSITY  OF  PARTICAL  IN  WHOLE  BLOOD 
C  SIGMSV  SCATTERING  CROSS  SECTION  OF  VEIN 
C  SIOMAV  ABSORPTION  CROSS  SECTION  OF  VEIN 
C  SIGMAS  SCATTERING  CRCXSS  SECTION  OF  TISSUE 
C  SIOMAA  ABSORPT  ION  CROSS  SECTION  OF  TISSUE 
C  OOEFS  SCATTERING  COEFFiaENT  TISSUE 
C  COEFSV  SCATTERING  COEFFiaENT  VEIN 
C  COBFA  ABSORPTION  COEFnClENT  TISSUE 
C  COEFA  V  ABSORPTION  COEFFICIENT  VEIN 
C  THETAC  CRITICAL  ANGLE  CALCULATE  FOR  SURFACE 
C  TISSUE  BOUNDARY 

CNO  INDEX  OF  REFRACTION  OUTSIDE  TISSUE 
CNT  INDEX  OF  REFRACTION  OP  TISSUE 
C  XO,  YO  COORDINATES  OF  CENTER  OF  CIRCULAR  VEIN 
C  RADIUS  RADIUS  OF  CIRCLE 
C  SLNGTH  ACCUMULATIVE  PATH  LENGTH 
C  NPHOTN  NUMBER  PHOTONS  TO  BE  INIECTED  INTO  TISSUE 
C INTIIALIZATION  OF  PARAMETERS 
DO  3001  -  1.100 
DO  300  J  =  1,100 
FLUX{J,I)=  0.0 

300  CONTINUE 
DO  30]  I  1,11 
DO301J  •  1,600 
ALR(U)  -  0.0 
ALT(I/)  -  OX) 

301  CONTINUE 

N-2 

AO)  -  0.0 

A(2)  =  0.0 

THETA  =  0.0 

PLNGTH  -  0.0 

RHO  -  5.0E03 

CALLSECNDS(ISBED) 

NPHOTN  =  )0000 
SIGMSV  *  56.58 
SIGMAV  -  0.542 
HEMAT  =  0.41 


ELSE 


CALL  TISSUE(RANDMl,RAm>M2,COEFA,COEFS,X,YJ»LNGTH, WEIGHT, THETA) 
ENDJF 

SLNGTH  «=  SLNGTH  +  PLNGTH 
IX  =  INT(AINT(X/10.0)) 
lY  =  INT(ABS(ArNT(Y/10.0)))  +  1 
ENOIP 

C  TERMINATION  OF  PHOTON  IN  TISSUE 

IF  (WEIGHT  IT.  XX»5)THEN 
WEIGHT  -  1.0 
GOTO  10 
BNDIF 

IF(y.LT.0D)THEN 

II^X  .GT.  49  .OR.  a  .LE.  -50  .OR.  lY  .GT.  100)THEN 

INC  -  INC  +  1 

ELSE 

FLUX(IX+514Y+1)  =  FLUX(IX+5UY+1)  +  WEIGHT 

INC*  INC+1 

ENDIP 

OOTPOl 

BL^ 

C  IF  PHOTON  IS  AT  SURFACE  INTERFACE  PROGRAM  CALLS 
C  SUBROUTINE  REFLECT  WHICH  DETERMINES  HOW  THE  PHOTON 
C  WILL  BE  REFLECTED 

CAlLREFLECr(N04n’,THETA,THETAC,Y4C,WT,WR,THETAI,THETAT 
t,RreRP41PAR4CS) 
lY  -  1NT(ABS(AINT{Y)))  +  1 
IXS  -  INT(AINT(XS/10j0)) 

IF(ILNGTH  JLT.  600)THEN 
ILNGTH*  INT(SLNGTH/10.0) 

ALT(14LNGTH)  -  ALT(1, ILNGTH)  +  WEIGHT*WT 
ALR(IJILNGTH)  -  ALR(1, ILNGTH)  +  WEIGHTWR 
lA  -  ABS(INTpCS/50.0))  +  2 
IFOA  .LB.  11)THEN 

ALT^OLNOTH)  -  ALT(IA,ILNGTH)  +  WEIGHT*WT 
ALR(U, ILNGTH)  =  ALR(IA,ILNQTH)  +  WBIGHT*WR 
ENDIF 
ENDIF 

C  RECORDING  OF  WHERE  PHOTON  HAS  TRAVELED  DURING  THIS  INCREMENT 

IF(IXS  .GT.  49  .OR.  DCS  .LT.  -SO)aOTO  111 

IFpX  .GT.  40  .OR.  IX  .LE.  -50  .OR.  lY  .GT.  100)GOTO  111 

IF(WT.EO.LO)THEN 

FUJX(1XS+5J,1)  -  FLUX(IXS+51,1)  +  WEIGHT 

GOTOlll 

ENDIF 

1F(WR.NE.J.0)THEN 

FLUX(IXS+5L1)  »  FLUX(IXS+5l.l)  +  WEIGHT'WT 
FLUX(IX  +51, lY  +1)  -  FLUX(1X  +51,IY  +1)  +  WEIOHT*WR 
ELSE 

FLUX(IX+514Y  +1)  *  FLUX(IX  +514y  +3)  +  WEIGHT 
ENDIF 

111  CONTINUE 


ENDIF 

GOTOl 

10  CONTINUE 

C  CONSTRUCnON  OF  OUTPUT  FILE  THAT  BUILDS  A  PILE  WHICH 
C  CONTAINS  INFORMATION  ON  THE  FLUX  DENSITY 

DO301Y  ■=  UOO 

WR!TB(1,’{"IV’.T4,/,»1Y  IX  FLUX’T)aY-l)*lO 
DO  20 IX  -  UOO 

WRITE(1.’(14,1X^4,1X,G13.6)’)(IY-1)*10,(K-51)‘J0.FLUX(IX.1Y)/ 

INPHOTN 
20  CONTINUE 

WRlTECl.’fEOF’)’) 

30  CONTINUE 

C  CONSTRUCTION  OF  FILE  THAT  OUTPUTS  INFORMATION  OBTAINED  FOR 
C  HISTOGRAMS  OP  TOTAL  PHOTON  PATHLENGTH  BEFORE  READMITTANCE 
C  AT  THE  TISSUE  SURFACE  INTERFACE 

WRITE(2,’(/,”OUTSIDE”./."LENGTH  ALL  D050  D50100”, 
nX,"D100150  D150200”)’) 

DO  303J»  1,600 

WRlTE(2,‘(6(G13Ji.lX))’)REAL{J)*10.0ALT(U)/NPHOTN 

!,ALT(2,1)/NPH0TN,ALT(3J)/NPH0TN.ALT(4,1)/NPH0TN,ALT(54)/NPH0TN 

303  CONTINUE 
WRITB(2.'0'*EOF')’) 

WRlTE{2,’("OUTSIDE2", /."LENGTH  0200250  D250300  D300350”, 

!1X,'T>3S04(I0  0400450  D450S00”)’) 

IM304K  •  1600 

WRITE(2.’(7(G13.6,lX))’)REAL(K)*10.0ALT{<y4)/NPHOTN, 

!ALT(7.K)/NPH0TN,ALT(8,K)/NPH0TN 

IAtT(»JK)/NPHOTN.ALT(n,K)/NPHOTN.ALT(ll,K)/NPHOTN 

304  CONTINUE 
WRITB(2,*fEOF’)’) 

WRrrE(2.’(/.”INSIDE”./.”LENOTH  ALL  D050  050100", 
I1X,"D100150D150200")’) 

DO  3051  -  1,600 

WRITE(2,'(6(ai3.6,lX))’)REAL(I)«10DAl-R(l.I)/NPHOTN, 

IALR(2J)/^PHOTN,ALR(34)/NPHOTN,ALR(4,I)/NPHOTNALR(5,I)/NPHOTN 

305  CONTINUE 
WRn'E(2.’("*EOF’)’) 

WR1TE(2,’(”INSIDE2", /."LENGTH  D200250  D250300  D300350”, 

II  JC,’’D330400  0400450  0450500”)’) 

D0306I«L600 

WRITE(2,’(7(013.6,]X))*)REAL(1)*10.0,ALR(0,I)/NPHOTN, 

!ALR(7,1)/NPH0TNALR(8,I)/NPH0TN 

I,ALR<9,?)/NPHOTN,ALR(10,I)/NPHOTN,ALR(I1,I)/NPHOTN 

306  CONTINUE 
WR1TE(2,’C'*E0F’)’) 

STOP 

END 

(.*»*«**.«.*....***o* . . . . 

(S  .*****.*«**.•«•**•.  reflection  . . 

c 

SUBROUTINE  REFLHCT(N0, NT, THETA, THETAC,Y,X.WT.WR, 
rrHETAl,THETAT,RPBRPJlPAR,XS) 

REALN0,NT 

THETAI  «  ABS{THETA  -  (3.1416/2.0)) 


CJ  O  r 


c 

c 

c 

c 

c 

c 
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TANTA  -  TAN(THETA) 

JPfTANTA  0.0)TANTA  =  0.00001 
XS-X-Y*(1D/TANTA) 

Y  --14)*Y 


dcgreei  >=  theUi*  180.0/(3.1416) 
degreet  <=  theuM80.0/(3.1416) 
degrcec  =  lbcttc*I80.0/(3.1456) 
J^lc(6/("«h««.”4a3.6/Hhrta=”^13-6/'^ 

f  lhet*c-",gi3,(^"i[o  >>^13.6  ”^i3.6)>) 

l<lflgreei,dcgrcet,y,degrMCA,xa 
*P(THerAJ  .LE.  THETAC)THBN 
™CTAT  »  AS1N((N0/NT)*S1N(THETAI)) 
1F(THETAI  .EQ.  0.0)THEN 
WT  -  1.0 
WR  =  0.0 
GO'TOIO 
ENDIP 


(M^ERP*RPERP  +  RPAR'RPAR) 


RPBRP 

RPAR  _ _ 

WR  -  (a5)*(RPERP*RPERP* 
WT-  LO.WR 


ELSE 

WT-ao 

WR-ID 


ENDIP 

return 

END 


REALNOJfT 

THETAI  -  ABSfTHETA .  (3.1416/2.0)) 

Y  ■  -TO'Y 

IPCTHETAI  IE.  THETAOl'HEN 
TOETAT  -  AS1N((N0/NT)*S1N(THETAI)) 

WR  .  {0.5)«(RPERP*RPERP  +  RPAR*RPAR) 
WT  -1.0.WR  ’ 


cc 
c 
c 
c 
c 
e 
c 
c 
c 
c 
c 

C  ELSE 
c  WT  -  0.0 

c  WR  =  IX) 

c  ENDIP 
c  RETURN 
c  END 

. 

VEIN  . . 

. . . . . . . . . . . 

TWOPI  -  2.0*(3.1416) 
theta  -  TWOPI*(RANDMl) 

^GTH  <*  (•li)*Al-OG(RANDM2))/(COEFSV) 

X  -  X  +  (PLNGTH'COS(TOETA)) 

Y  ..  Y  +  (PLNGTH».VJN(THETA)) 

R^RN  " 

END 

. . . . 

tissue  . . . 
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SIGMAS  >  56^ 

S1GMAA  -  0.0001 

COEFS  -  (RHO*SIOMAS*HEMAT*(1.0'HEMAT))/(1.0E06) 

COEFSV  -  (RHO»S10MSV*HEMAT*(1D-HEMAT))/(1.0E06) 

TWOPI  =2D*(S.1416) 

COEFA  -  (RHO*SIGMAA)/(1.0E06) 

COEFAV  «  (RHO*SIGMAV)/(1DE06) 

W0GHT  -  ID 
WR  -0.0 
WT  -OD 
NO  -  1.0 
NT  -155 

THETAC  ■=  ASIN(N0/VT) 

RPAR  -  0.0 
RPERP  =  0.0 
XO  -OD 
YO  -  -200D 
RADIUS  -  150.0 

C  INITIALIZATION  OF  RANDOM  NUMBER  GENERATOR  AND  THE 
C  GENERATION  OF  PHOTONS  TO  BE  INSERTED 

1  CALLHSRPST(1SEED) 

DOIOI-  LNPHOTN 
X  ~OJ0 

Y  -OD 
IX  -0 
lY  -0 
INC-  1 
SLNGTH  -  OD 
SWITCH  -  1.0 

FLUX(S1,1)  -  FUJX(51,1)  +  10 

C  CVOXXJLATION  OF  RANDOM  NUMBER  AND  POINT  WHERE  ITERATION 
C  OF  PHOTON  PATHLENGTHS  WHILE  IN  MEDIA 

CALLHSRPON(N^) 

IF  (INC  .GT.  5000)  GOTO  10 
RANDMl  -  A(]) 

RANDM2  -  A(2) 

C  INSERTION  INTO  TISSUE 

IF  (SWITCH  .EQ.  1.0)THEN 

x-ao 

PLNGTH  -  (-1.0*AIjOG(RANDM2))/(CXJEFS) 

Y  -  Y  +  (.1.0)*PLNGTH 

swrrcn-ao 

ELSE 

C  CALCULATION  OF  PHOTON  POSITION  RELATIVE  TO  VEIN  AND  TEST 

C  TO  SEE  WHETHER  PHOTON  IS  INSIDE  OR  OUTSIDE  THE  VEIN,  PHOTON 
C  ENTERS  APPROPRIATE  SUBROUTINE  WHERE  PATHLENGTH  AND  ABSORPTION 
CIS  CALCULATED 

RORC  -  SORT  ((X  -  XO)  *  (X  -  XO)  +  (Y  -  YO)  *  (Y  -  YO)) 

IF(BaRC  .LE.  RAD1US)THEN 

CAU.VE1N(RANDM1JIANDM2,C0EFAV,C0EFSV,X,Y4>LNGTH, WEIGHT, THETA) 


noon*' 


\r‘\ 

. . . . 

SUBROUTINE  T1SSUE(RANDM1JIANDM2,C0EFA,C0EPSAY.FLNGTH 
!, WEIGHT, THETA) 

TWOPl  =  2j0*(3.1416) 

THETA  -  TWOW*(RANDMl) 

PLNGTH  =  (-LO»ALOri(RANDM2))/(COEFS) 

X  -  X  +  (PLNGTH*COS(THETA)) 

Y  -  Y  +  (PLNGTH*SIN(THETA)) 

WEIGHT  »  WEIGHT»EXP((-1.0)*(PLNGTH)*(COEFA)) 

RETURN 

END 


sccnds 


SUBROUTINE  SECNE>S(JSEC) 
inlcgcr*2  idat(6) 
CSinclude;7sys/ins/timejiu.(tn’ 

C  call  CAL  SDECODB  LOCAL  TlME(ldat) 

C  lSEC-lNT(idat(6)  +  10*idat(^) 

C  return 
C  end 


